NUMERICAL INVERSE SCATTERING FOR THE TODA LATTICE 


DENIZ BILMAN AND THOMAS TROGDON 


Abstract. We present a method to compute the inverse scattering transform (1ST) for the famed Toda lattice 
by solving the associated Riemann- Hilbert (RH) problem numerically. Deformations for the RH problem are 
incorporated so that the 1ST can be evaluated in 0(1) operations for arbitrary points in the (n, t)-domain, 
including short- and long-time regimes. No time-stepping is required to compute the solution because (n, t) 
appear as parameters in the associated RH problem. The solution of the Toda lattice is computed in long-time 
asymptotic regions where the asymptotics are not known rigorously. 




o 

(N 

Oh 

<D 

C/3 

oo 

<N 


C/3 

d 


> 

OO 

oo 

r- 


oo 

o 

in 


X 


1. Introduction 

We consider the numerical solution of the Cauchy initial value problem for the doubly-infinite Toda lattice 


( 1 ) 


d t a n {t ) = a n (t)(b n+1 (t) - b n (t )) 
dtb n (t) = 2 (a n (f ) * 2 - a n _i(f) 2 ) 


a n ( 0 ) = a° > 0 , and 6 n ( 0 ) = 6 °, 

for (n, i) G Z x M with solution^ (a(f), 6(f)) G ^(Z) x £(Z) satisfying 
(2) £ o(n) (|a n (t) - 5 1 + \b n {t)\) < oo, for all f G ] 


a(n) = e' 5|n| , 


nEZ 


for somc 0 S > 0 . 



2000 


0.65 



0.55 


2300 


Figure 1. An example solution of the Toda lattice computed at t = 2000 with the method 
presented here. The initial data is given by a n (0) = 11/2 — ne~ n +Tl | and b n ( 0) = nsechn. This 
initial data produces dispersive radiation (center panel) and four solitons, two traveling each 
direction (left and right panels). 
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"*We omit subscripts to refer to the functions defined on Z. 

2 Many results for the Toda lattice hold with less restrictive choices for a(n) [31]. For example, the inverse scattering transform 
method described below can be applied for data in the so-called Marchenko class (i.e., cr(n) = 1 + |n|). We impose exponential 
decay for the convenience of the numerical implementation. 
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The Toda lattice was introduced by Morikazu Toda in [36] (see also [35]). The Toda lattice is a completely 
integrable model for a one-dimensional crystal. The form ([Tj) we use in this paper is the Toda lattice written 
in Flaschka’s variables [15] (see also the work of S. V. Manakov [22] 1. This system has been studied in great 
detail because it is the prototypical discrete-space, continuous-time infinite-dimensional integrable system. 

A consequence of the complete integrability of the Toda lattice is an associated inverse scattering transform 
method (ISTM). The ISTM first maps the initial data to a spectral plane where its time evolution is simple 
via a transformation called direct scattering , see Section 2.2 Then, at a given time t, the evolved spectral 


data is mapped back to the physical plane to find the solution values a n (t ) and b n (t) by a transformation 
called inverse scattering. This inverse problem is solved by considering an associated oscillatory Riemann- 
Hilbert (RH) problem. RH problems are boundary-value problems in the complex plane for sectionally 
analytic functions. General references are [2] El El E31 EH So] ST]. From an analytical point of view, the 
benefit of studying the RH problem is that asymptotics can be extracted by the method of nonlinear steepest 
descent. Broadly, the method works by deforming the contours of the RH problem as in the classical scalar 
method of steepest descent to turn oscillatory terms to exponentially decaying terms. See mmm for some 
implementations of this method. 

As for specific applications of the method of nonlinear steepest descent to the Toda lattice, we refer 
the reader to the work of Kamvissis m and Kruger and Teschl PSEU. These works give explicit long¬ 
time asymptotics for the solution of the Toda lattice in the soliton and dispersive regions that we define 
in Section 2.4 The work of Kamvissis gives the asymptotic behavior in the Painleve region (as defined in 


Section 2.4) for non-generic initial data. 

We approach the ISTM from a numerical perspective. See Figure [l] for a sample solution computed with 
our method. Given initial data with sufficient decay (see ([2])) we are able to compute the solution at a given 
point (n, t), to a given accuracy, in a bounded number of operations by solving the RH problem numerically 
and incorporating the deformations used in the method of nonlinear steepest descent. Stated another way, 
for any given (n, t) we give an 0(1 ) algorithm to compute the two values {a n (t), b n (t)}. No time-stepping 
is required to compute these values. This methodology has been previously applied to the KdV and mKdV 
equations [42], the focusing and defocusing NLS equations [4.3] . the Painleve II equation j28[:30j and orthogonal 
polynomials on the line (30j. We compute the solution of the Toda lattice for arbitrarily large values of t. The 
complexity and the accuracy of the methodology is discussed in m- The code we have developed is available 
at [44] , With the current state of the art, the use of deformations appears to be necessary as the associated 
RH problem is increasingly oscillatory as \n\ or t increase. Deformations can be avoided in some cases with 
oscillatory integral techniques [38] but unfortunately these methods are not currently general enough to use 
in the case of the Toda lattice. 

The RH problem associated with the Toda lattice has fundamental differences from problems previously 
solved numerically. First, the fundamental domain is the unit circle as opposed to the real axis as in the cases 
of the KdV, mKdV and NLS equations. Second, the RH problems in the previously solved cases have had their 
deformations worked out in detail. In this paper, we develop deformations for the Toda lattice in regions of 
the (n, t)-plane where deformations do not exist in the literature including the determination of the so-called 
^-function, see Appendix [C] We believe our deformations will be required for the future asymptotic analysis 
of the Toda lattice. Importantly, the asymptotic regions we define, and the deformations performed therein, 
cover the entire (n, t ) plane. This is something, to our knowledge, that has not been performed previously 
in the literature. Finally, we encounter some interesting technical challenges in computing the functions used 
in the deformations, see Appendices [A] and C.3 Also, in light of the work in [5] we believe this numerical 
method will be useful in studying non-integrable, Hamiltonian perturbations of the Toda lattice. 

In this paper we do not consider large amplitude data. Large amplitude data induces singular behavior 
in the solution of the Toda lattice akin to the behavior in the small dispersion limit for the KdV equation, 
for example. This affects numerical methods in a critical way. Although, we do not consider large amplitude 
data, throughout the manuscript we include footnotes that highlight the complications that arise for larger 
initial data. We also do not treat the case where poles in the Riemann-Hilbert problem are very close to the 
unit circle. The methodology described here can be used to handle this case accurately with some extra work 
and deformations. 

The paper is organized as follows. In Section [2] we give the background material on the ISTM for the 
Toda lattice. The direct scattering and inverse scattering maps are discussed along with a discussion of the 
(asymptotic) regions (soliton, dispersive, Painleve, collisionless shock, and transition) of the Toda lattice. We 
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also describe the fundamental deformations, which are performed in all of these regions. Section [3] provides 
a step by step guide outlining the computational procedure for obtaining the solution of the Toda lattice. 
The majority of the paper is devoted to Section [4] where we discuss, and explicitly derive, the deformations 
of the RH problem in each region. In Section [5] we discuss the numerical solution of RH problems. Finally, in 
Section[ 6 ]we give some numerical results including an error analysis. We include five appendices. Appendix |A| 
discusses the numerical solution of singular, but diagonal, RH problems. Appendix [B] gives a deeper discussion 
of computing the eigenvalues of Jacobi operators. Appendix [C] details the ^-function that is used in both the 
collisionless shock and transition regions. Appendix [D] contains the vanishing lemma and a discussion of the 
unique solvability of RH problems considered in this work. Lastly, Appendix [E] gives a proof that Jacobi 
matrices whose reflection coefficient attains the value —1 at the edges of its continuous spectrum forms an 
open dense subset of the Marchenko class (c.f., <r(n) = 1 + |n| in Q) of Jacobi matrices. This implies that 
for an open dense set of initial data, the long-time behavior of the solution of the Toda lattice exhibits a 
collisionless shock region: see Section 2.4 below (see also ISIS!)- If the reflection coefficient does not attain 
the value —1 at the edge of the continuous spectrum then the collisionless shock region is absenl|^| 


2. Background material 

We use this section to cover theoretical background and fix notation. 

2.1. Integrability and Lax pairs. The complete integrability of the Toda lattice was proved by H. Flaschka 
in 1974 in a sequence of papers m and ESI, and independently by S. V. Manakov in |22j . Introduce the 
second-order linear difference operators L and P defined on f 2 (Z) by 

(3) (-L/)n — Q"n—lfn—l T b n f n T 0-nfn- 1-1, 

(4) (F > /) n — CL n —if n —i T CL n f n -|_i . 

and note that in the standard basis L = L({a n } ne z, {frnjnez) is a Jacobi matrix (symmetric, tridiagonal with 
positive off-diagonal entries) and P is a skew-symmetric matrix, i.e., P T = —P: 




bn—1 C^n— 1 

0 

\ 



0 d n —\ 

0 

\ 

L = 


tin— 1 bn 

tin 


and P = 


^n—1 0 

On 




0 d n 

bn+ 1 




0 d n 

0 



V 



‘•J 


V 





The system of equations given in ([!]) is equivalent to 

d t L(t) = [.P(t),L(t )] = P(t)L(t) — L(t)P(t), 

and (P. L ) is called a Lax pair. Its existence shows the complete integrability of the Toda lattice. A conse¬ 
quence of complete integrability (or of the Lax pair) is the existence of an inverse scattering transform for the 
Toda lattice. 


2.2. Direct scattering: definition of the scattering data. Since L is a bounded self-adjoint operator 
the spectrum <r(L) C M. Furthermore, ([ 2 ]) implies that the spectrum of L consists of a purely absolutely 
continuous (a.c.) part 

Cac(A) = [-1, 1], 

and a finite simple pure point part 


Opp(L) = {Aj : j = 1, 2,..., N} C (-oo, -1) U (1, +oo). 
For convenience we map the spectrum via the Joukowski transformation: 

A = \ { z + z 1 ) , ^ = A — \j A 2 — 1, A 6 C, \z\ < 1. 


-This also implies that the transition region, defined in Section 


2.4 


is also absent. 
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Here the square root \/A 2 — 1 is defined to be positive for A > 1 and cr ac (-L) is the branch cut. Under this 
transformation, the a.c.-spectrum, [—1,1], is mapped to the unit circle T and the eigenvalues A j are mapped 
to Cf\ with (j € (— 1 , 0 ) U ( 0 , 1 ) via 

( 4 5 ) = \ (Cj + C ) , 

for j = 1,2,..., iV. For any z with 0 < \z\ < 1, z ±1 the equation 

( 6 ) Up = 

has two unique solutions, <p + and </?_, normalized such that 


(7) lim z Tn (p±(z; n) = 1. 

n—>-=boo 

For fixed n, <p±{z\ n ) are analytic functions of z, 0 < \z\ < 1. With the assumption Q of exponential decay in 
the initial data, the functions <p±{z\ n) extend analytically to 0 < \z\ < l + r(<5) where r(8) > 0 depends on the 
decay rate of a(n) in ([2]). It follows from Green’s formula that the Wronskians W ( <p±{z ; ■), ! p±(z~ 1 ', •)) are 
independent of n, and evaluating them at ±oo we observe that {(p±(z\ ■),'p±{z~ 1 \ •)} are two sets of linearly 
independent solutions. We define the transmission coefficient T(z) and the reflection coefficients R±{z) by 
the scattering relations for \z\ = 1 

^ T(z)<p+(z,n ) = R-(z)<p-(z;n) + p _ (.z _ 1 ;n) , 

T(z)ip-(z, n) = R+(z)ip + (z-,n) + ip+ (z^;n) . 

Also, for general data the transmission coefficient has a meromorphic extension inside the unit disk \z\ < 1, 
with finitely many simple poles at Cj, |Cj| < 1, j = 1,2 ,..., N. The residues of T(z ) are given by: 

(9) Res T(z) = -Cjl+jUl 1 = -Cjl-jVj, 

Z=Cj J 


where 

. x 1 

( 10 ) 7 ± j = - k - 

J lb±(0; ■)&(*) 

are the norming constants and pj is the associated proportionality constant: tp-(Q, •) = /J,jip+(Cj, ')• Due to 
the assumption on the initial data ([ 2 ]), the relations Q remain valid in an annulus containing the unit circle 
and therefore R±(z) and T(z) are meromorphic in this annulus. 

One reflection coefficient, one set of norming constants, and the set of eigenvalues is sufficient for recon¬ 
structing L via the inverse scattering transform for Jacobi matrices whose coefficients decay sufficiently fast 

[31] • Define 

R{z) = R+(z) and jj = j +jj , j = 1,2,... ,1V, 

and the set 

S(L) = {i?(z),{0}f =1 ,{7,}f =1 } 

to be the scattering data for the Lax operator L. For a more detailed account of the scattering theory for 
Jacobi matrices, see [3l] or [21], 


2.3. Inverse scattering: the Riemann—Hilbert problem. We phrase the inverse problem in terms of a 
sectionally meromorphic RH problem. In what follows, plus (+) and minus (—) sides of a contour correspond 
to the left and right sides by orientation, respectively. And m±(z)= m±(z) denote the boundary values of a 
function m(z) as z tends to the relevant contour from the ± side. 

RH Problem 1 . Let the unit circle T have counterclockwise orientation. As in m, we seek a function 
m: C\T -A C lx2 that is sectionally meromorphic, continuou^j up to T, with simple poles at Cf 1 , 3 = 1 • • •, IV, 
and satisfies: 


4 Throughout this paper, unless we specify otherwise, we look for solutions of the Riemann-Hilbert problems that are continuous 

up to their jump contours. 
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the jump condition: 

m + (z;n,t) = m~(z;n,t)J(z;n,t), z£ T, J(z;n,t) = 


0 


l-\R(z)\ 2 -R(z)e~ e ^ n ^ 

R(z)e 9 ^ n ’ t '> 1 


(ID 


• the residue conditions: 



Res mlz; n, t) 

= lim 777(2; n 

;t) ( 


2_> C j 

Res m(z:n,t ) 

= lim 777(2; 

77, t) 


^c- 1 

• the symmetry condition: 




777 (2 _1 ; 77, 

*) = 




0 C7 1 7je 0( ^ ;n - t) 


'0 1 

.1 0 




• the normalization condition: 

(12) lim m(z; n, t) = (mi m 2 ) , with m\ ■ m 2 = 1 and m± > 0 . 

Z—>00 V J 

Here the exponent 9(z; n, t ) in the jump matrix J is given by: 

(13) 9(z\ n, t) = t (2 — 2 _1 ) + 2nlog(-sr), 

The symmetry condition ensures that RH Problem [l] has 


and it can be shown that R(z) = R(z x ) 
unique solution for all values of (n,t) (see Section 3 in [20]). 

Remark 2.1. The associated matrix RH problem for RH Problem [I] (the RH problem with a 2 x 2 unknown 
function satisfying the same jump condition, normalized to the identity matrix at infinity and no symmetry 
condition, see Definition 5.1 below) may not have a solution for some exceptional (n,t) values. Indeed, these 
exceptional values are guaranteed to exist when N 0, see m Lemma 2.6] and the preceding discussion, for 
example. But as illustrated by this example, such an exceptional value of (n, t) occurs when ~ 

|("l — C 1 ], i.e. near the peak of a soliton. This phenomenon is a consideration in the numerical method 


developed in this work, see Remark 6.1 


We have the following well-known and important fact: 

Proposition 2.2. For generic initial data R(±l) = —1 and hence m|(—1) = mij (—1) = 0. If the potentials 
(a 0 — 1/2,6°) tend to zero exponentially as \n\ -A 00 , R(z) is analytic in a neighborhood of T. Moreover, 
m ± (z) have analytic extensions across T and hence mf(z) and mf(z) have a zero of at least first order at 
z = - 1. 


Proof. First, that generically ii(±l) = —1 is shown in Appendix |E} Let rn(z) be the solution of RH Problem[lj 
Then at z = — 1 

m + (-l) = m“(—1) 

because R(— 1) = — 1. The first component of this equation gives 

1 ) = -mf(-l). 

But the symmetry condition ( |11[ ) gives that mf(— 1) = mf (—1) so that m/"(—1) = mf (—1) = 0. Analyticity 
follows from considering the Volterra summation equations ( |92| ) which forces the zero to be of at least first 
order. □ 


We proceed with a lemma for recovering the potential from the unique solution of RH Problem [TJ 

Lemma 2.3. The solution ( a(t),b(t )) E £(Z) x £(Z) to the initial value problem (jTj) for the Toda lattice can 
be recovered from the asymptotic behavior of m(z\ n, t ) near z = 0, 

m(z] n, t ) = ^A n (t) (l — 2B n -\(t)z) ^(t) + 2 B n (t)z)^ + O (z 2 ) , 
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t 



or z = oo, 

(14) m(z; n, t) = (a^)(1 + 2B n {t)z~ l ) A n (t)(l - 2B n _ 1 (t)z~ 1 )^J + O (z~ 2 ) 

where A and B are defined by 

OO OO 

A n (t) = 2 a n (t) and B n (t) = - ^ bj(t). 

j=n j=n+l 


Furthermore, (12) ensures A n (t) > 0 and hence a n (t) > 0 for all t > 0. 


2.4. Asymptotic regions. In this section we discuss asymptotic regions for the long-time asymptotics of 
the Toda lattice with decaying initial data. A rigorous study of long-time asymptotics for solutions of the 
Toda lattice equations was recently carried out in [20] and [21] in the soliton and the dispersive regions (see 
below), but the question of long-time asymptotics in the region |n|/t ~ 1 has not been addressed in generality 
so far. The long-time behavior of solutions in this region was studied in m under the additional assumptions 
that no solitons are present, that is, the RH problem has no poles, and that |i?(±l)| < 1. Under the latter 
assumption, the solution is given asymptotically in terms of a Painleve II transcendent in the region \n\/t ~ 1, 
[19]. However, one generically has R(z = ±1) = —1 (We give a proof of this fact in Appendix [E]). In this 
case, an additional region called the collisionless shock region appears as the stationary phase points of the 
jump matrix coalesce at z = ±1, and one needs to introduce additional contour deformations, employing the 
so-called ^-function method, to bridge the dispersive and the Painleve regions. In the current work, we present 
new deformations of the associated RH problem for this unstudied region \n\/t ~ 1 with generic initial data. 
These deformations are essential to compute solutions numerically. 

Introduce constants, Cj > 0, to divide asymptotic regions. 

1. The dispersive region. This region is defined for |n| < c±t, with 0 < ci < 1. Asymptotics in this region 
were obtained in [211] . 

2. The collisionless shock region. This region, to the best of our knowledge, has not been addressed in 

the literature. It is defined by the relation c\t < |n| < t — (logf) 2//3 . Asymptotics are not known 

in this region. 

3. The transition region. This region, to the best of our knowledge, is also not present in the literature. 
The region is defined by the relation t — cfi^^ 3 (log t) 2//3 < |n| < t — C3f 1 / 3 . Asymptotics are not known 
in this region. An analogue of this region was first introduced for KdV in [42]. 

4. The Painleve region. This region is defined for t — cfi 1 ^ < |n| < t + cfi 1 ^. Asymptotics in this region 
were obtained in |19| in absence of solitons and under the additional assumption that |-R(z)| < 1. 

5. The soliton region. This region is defined for |n| > t + cfi 1 ^. Let > 1 denote the velocity of the 
k th soliton and choose u > 0 so that the intervals (Vk — v,Vk + v), k = 1, 2,..., N, are disjoint. If 


\ n /t~Vk\ < v, the asymptotics in this region were obtained in |20j and [21] , It will follow in Section 4.2 
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that the deformation in the Painleve region for |n| > t is the same as that for soliton region although 
one should expect the long-time behavior to be different in each region. We will see that from a 
numerical perspective the Painleve region for \n\ > t can be identified with the soliton region. 

2.5. Fundamental deformations of the inverse problem. As outlined in Section [3] below, the procedure 
for numerical solution of RH Problem [I] involves a sequence of deformations, dictated by the values of ( n,t ), 
which result in a numerically tractable RH Problem satisfied by a sectionally analytic vector-valued function. 
In this section we present two deformations that are performed for all values of the parameters ( n,t ). Our 
first step is to remove the poles (if any) from the sectionally meromorphic RH Problem [I] This is achieved 
by introducing small circles centered at each pole and using the appropriate jump conditions on these new 
contours [8] . Fix e > 0 such that 


(15) 


e < — min < min {|Cj|},min {||Cj| - 1|} , min {|Cj - Cz|} 


1 <j<N 


1 <j<N 


1 <j,l<N 

i¥* 


and define the circles Df by 


Df = {*: 


.±1 


-Ci|=e}> 3 = 1,2,..., AT. 

± 


The choice (pL5|) of e guarantees that the disks enclosed by the circles D~j, j = 1,2, 


, IV, do not intersect 


each other or the unit circle and none of them contains the origin. We define fh{z\ n, t ) by 


(16) 


m(z; n, t) := < 



N, 


z ~ Cjl < e, 3 = 1 , 2 ,.. 

/T 1 - Q | < e, j = 1,2,..., N, 


^0 1 

m(z;n,t), otherwise. 

It is straightforward to show that rh(z) solves the following sectionally analytic RH problen^j 

RH Problem 2. 


fh + (z; n,t) = < 


(17) 


m ( z',n,t)J(z',n,t ), 

( 1 

m (z] n, t ) 


fh ( z\n,t ) 


z-Cj 

6(Cn \r. 

I z Zj e J 


z-C~ 


.0 


0 1 
1 0 


1 

\z\ > 1 , 


zGl, 


21 G Df, j = 1 , 2 ,. 

..,N, 

2 e Dj, j = 1 , 2 ,. 

..,N, 


fh (z 1 ; n, t) = fh(z; n, t) 
lim fh(z\ n, t ) = (mi m 2) , with mi ■ m 2 = 1 and mi > 0 , 

Z —¥OO V 7 

where {Dj} are oriented counter-clockwise and {Dj} are oriented clockwise. 

This deformation brings in a possibility of exponential growth for t > 0 in the new jump matrices on {D±}. 
There are two cases to distinguish. If R eO(Cj',n,t) < 0 the jump matrices introduced in RH Problem [2] have 
exponential decay to the identity as t — > 00, which is what we desire. If Re 0(Q] n,t) > 0 for some (j, however, 
the jumps around such poles are unbounded as t —> 00. Following the approach in [8] (see also [21]) we employ 
a conjugation procedure to restate our problem so that the jump matrices tend to the identity exponentially 
fast as either |n| or t tend to infinity. Let K U) t C {1, 2,..., N} denote the index set for (j (if any) such that 
[7j(t)| exp (Ref?(£j; n, t)) > 1. We set 


^From here on we state RH problems only in terms of their jump condition, jump contour, symmetry condition, and 
normalization. 
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(18) 


Q(z) = 


n 


z ~Cj 

-i 


\ 


j£K n , t Z 

0 


II 


V i&K n , t z y 

and Q(z ) is naturally defined as the identity matrix if K n y = 0. For j £ K ll} t define 

2-C?_ 

:n.t) \ 

Q(z), 


(19) 


m(z ; n,t) = < 


1 


m{z\ n, t ) 


m(z; n, t ) 


C me 


6((j\n,t) 


JjCje 1 




z - (A < e, j e K , 


n,£ 5 


2-0 

0 

2 0 1 


0 

20' 

20- 1 


Q{z) 


zCjlje 


B(( 


The matrices 


1 


m(z-,n,t)Q(z), 

2-0 


12 0| < <0 j £ Kn,ti 

otherwise. 


0'7je' 




2-0 


-iO- 9(c ^' M) 

0 


Q(s;) and 


0 

2Q— 1 


z 0'7j e 


S(C j;n,t) 


zCj ' 

20— 1 

1 


<20) 


have removable singularities at Q and £• , respectively. Now, 'm(z) satisfies: 


RH Problem 3. 


m + (z; n, t) = < 


( 20 ) 


m ( z\n,t)Q 1 (z)J(z-,n,t)Q(z), 


m ( z\n,t)Q l (z) 0 : T? 


2 0 


,0 


1 


m ( z;n,t)Q L (z) gQ-i 


,20 7:/ e ( 


,®(C ,■;«,*) 


2 6 T, 

Q(yZ)i Z £ Dj , J £ Kn,ti 

I Q(^), 2 £ D~, j £ iO n ,t, 


m 


/ i 0\ 

( z;n,t)Q x (z) I J Q( z )> z e -Dj , j ^ ^n,t 


m ( z;n,t)Q *(2) 


'j <X,e 6( ^ n/) 


,0 


zC -»' -1 I Q(z), z e Dj , j £ K n . t} 


z—> oo 

where m,- £ M are given in 


m(0;n,t) = m(oo;n,t) ^ ^ Q(0). 
lim m(z;n,t) = (mi m2) , with mi • m2 = 1 and ?bi > 0 , 

>• - i.r'v'l ' ' 


Note that we have relaxed the global symmetry condition © (required to hold for all 2 £ C) present 


m 


and RH Problem [2] to an asymptotic symmetry condition (required to hold only at 2 = 00) in 
As we shall see in Section [5j doing this does n ot c ause a problem for numerical solution of 
and it simplifies some calculations, see Remark 5.7 If N is largqu the entries in these jump 


RH Problem 1 
RH Problem 3 
RH Problem ^ 

matrices can a so be large for finite (n, t) even though we have decay to the identity as \n\ or t becomes large. 
We consider small values of N. Finally, to simplify the notation, let Xj t ± denote the jump matrices that are 
used to invert the exponentials for j £ K n j. 

(^ __\ l 1 

( 21 ) 


Xj + (z;n,t) = 07 ,-e ^ 


0 (C,■;«.*) 


,0 


1 


and Xj-(z‘,n,t ) = zQ-i 




®One would expect large N if the data has large amplitude. 
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and Yj t ±(z;n,t) denote the jumps introduced in RH Problem [3] for j / K n p 


( 22 ) 


Y j,+(z;n,t ) = ^,- 7i e 9(c J ;n,t) -jJ and Yj_(z;n,t) = 


' Kive e( ^ n ’ tr 

*0-i 


,0 1 / 

To summarize, for all values of ( n,t ), we initially perform the following chain of deformations: 

[m(z\ n, t); RH Problem [lj i —> [fh(z ; n, i): RH Problem [2j i— > [ m(z ; n, t); RH Problem [3]. 


3. The Methodology: A step by step guide 


In order to compute solutions of the Toda lattice we must perform each of the procedures oulined in 
Section [2] numerically. Here we outline what this entails and give a brief discussion of each step. The sections 
that follow describe many of these steps in an increasing level of detail. We perform the following: 

3.1. numerical computation of the scattering data, 

3.2. deformation of the vector RH problem, and 

3.3. numerical solution of the deformed vector RH problem. 

Recall that we are computing with initial data that has exponential decay as described in ([2]). 


3.1. Numerical computation of the scattering data. 

• Computing R(z). 

For z in a neighborhood of T, on which R(z) is analytic, we look for solutions of Lip = ^ (z + z -2 ) <P 
which behave like z ±n as n A ±00. We define two new functions f + (z',n) = ip + {z\n)z~ n and 
f-(z;n) = p-(z',n)z n so that we have f±(z;n ) A 1 as n A ±00. Then f± satisfies 

a n -iz Tl f±{z\n - 1) + (b n - \ {z + z -1 )) f±(z\n) + a n z ±l f±(z\ n + 1) = 0. 


This can be effectively solved using back substitution on n > 0 using the appropriate boundary 
conditions for f+(z\ •): for K large approximate /+ by using the condition f+(z; K+l) = f+(z ; K) = 1. 
The constant K is chosen so that \b m \, \a m — 1/21 are both less than machine accuracy for \m\ > K. 
For n < 0 a similar method works for f-(z\ •) by setting f-(z\ —K — 1) = f+(z m , —K) = 1. Matching 
these approximate solutions at n = 0 yields an approximation of the reflection coefficient. 
Computing {CfS C^\ ■ ■ ■, C ^ 1 }• 

Computing is equivalent to computing the ^ 2 (Z) eigenvalues of the doubly-infinite Jacobi matrix L 
that is defined in ([3]). We approximate the eigenvalues of L by computing eigenvalues of a Lk that are 
outside the interval [—1,1] for a large value of K [5]. This method might fail to capture eigenvalues of 
L that are close to its continuous spectrum. We present an example illustrating this case and provide 
the underlying spectral theory for Jacobi matrices in Appendix [B] We check whether we successfully 
capture all of the eigenvalues by computing the inverse scattering transform at t = 0 and comparing 
the reconstructed solution to the initial data. If these values differ by a user-prescribed tolerance, we 
employ Newton iteration to compute the (simple) zeros of 1 /T(z), since, as mentioned in Section 2.2 


are the (simple) poles of the transmission coefficient, T(z). 


Computing { 71 , 72 , • • • ,7Ar}- 

At the points z = (j, the solutions ip±{z\ •) are proportional and thus both lie in £ 2 (Z) with exponential 
decay asn-> ± 00 . We compute the norming constants by computing n) for n > 0 and <p- ((j ; n) 

for n < 0. We determine the proportionality constant pj (see <§) by matching these solutions at n = 0. 


We then calculate the £ 2 -norm ||^j + (£j; Oll^z) piecewise, using ip + ((j-,n) for n > 0 and pj 1 ip-((j-,n 
for n < 0. Thus, by recalling (10), we obtain the norming constant 7 y. 

lj = 011 ^ 2 ( 2 ) • 


3.2. Deformation of the vector RH problem. For most values of (n, t), RH Problem[3]has high-oscillation 
in the jump matrix J(z; n, t). To be able to accurately compute the function m(z) (or m(z)), the RH problem 
needs to be deformed to control these oscillations. This is the objective of the Deift-Zhou method of nonlinear 
steepest descent (see, for example, m ei mu im). The input to deform the vector RH problem is both the 
parameter set (n, t) and the numerically computed scattering data. 
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Choose the (asymptotic) region for ( n,t ). 

For each paiin ( n,t ), n,t > 0, we need to associate a region. This region will dictate how to deform 
the RH problem. We use the five regions introduced in Section 2.4 and one additional region: 

0. a region where no deformation is made. 

In our code, we use tom [l] to choose the region: 


Algorithm [T) Choosing a region. 

Data: (n. t ) 

Result: The region for deformation. 

set c i = .96; set C2 = .2; set C3 = 3/2; set C\ = 3; set po = y/l — (n/f) 2 ; 
if n 2 + 4f 2 < C\ then 
no deformation 
else if n> t then 

use soliton region deformation 
else if t — C3 (f 1 / 3 + 1) < n < t + + 1) then 

use Painleve region deformation 
else if n/t < c then 

use dispersive region deformation 
else if n/t < 1 and — 2 log po/ftp^) < co then 
use collisonless shock region deformation 
else 

use the transition region deformation 

end 


Notes ON Algorithm 1. This is the algorithm for choosing a region. See (29) for the 
appearance of po in the analysis. The constants defined in the algorithm are user specified and 
can be adjusted on case-by-case basis. In our code we leave them fixed as displayed. 


• Deform the jump contours and compute auxilliary functions. 

Once the region has been chosen the deformation of the vector RH problem has to be implemented. 
This usually follows the ideas of the Deift-Zhou method of nonlinear steepest descent applied from a 
numerical perspective For each output of Algorithm [l] one has to determine and hard-code appro¬ 
priate deformations which often introduce additional, or auxilliary, functions that must be computed. 
These deformations and auxilliary functions are discussed in great detail in Section |4j 


3.3. Numerical solution of the deformed vector RH problem. Once the deformed vector RH problem 
is in hand, we can proceed with its numerical solution. We must perform the following steps: 


Compute the solution of the associated matrix RH problem. 

The vector RH Problem [3] and its deformations have a normalization at infinity that is difficult to treat 
numerically. So, we numerically solve the associated matrix RH problem whose solution is defined to 
be a 2 x 2 matrix-valued function that has the same jump condition, has no symmetry condition and 
is normalized to tend to the identity matrix at infinity. As discussed in Remark 2.1, this matrix RH 


problem may fail to have a solution at some exceptional values of ( n , t ) yet we can still reliably solve 


it numerically, see Remark 6.1 


Construct the solution of the deformed vector RH problem. 

Because the rows of the associated matrix RH problem are linearly independent the solution of the 
deformed vector RH problem must be a linear combination of the rows of the associated matrix RH 
problem. The combination is determined by solving a 2 x 2 eigenvalue problem. 

Extract the solution of the Toda lattice at (n,t). 

Once the solution of the deformed vector RH problem is computed the deformations must be reversed 
and a Taylor expansion is performed to compute A n (t ) and B n (t ) using Lemma 2.3 This procedure 


^We only consider t > 0 and n > 0. A transformation is used to treat n < 0, see Remark 


4.1 
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is repeated for (to + l,t). Using A n (t), A n+ i(t), B n (t) and R n _i(i) we have 


(23) 


a n (t) = -A n (t)/A n+ i(t), b n (t ) = B n (t) - B n _i(t). 
4. Deformation of the vector RH problem 


In this section we present the deformations of RH Problem [3] which are required in each asymptotic region. 
These deformations involve explicit functions that need not satisfy the global symmetry condition present in 
RH Problem [I] (or RH Problem [3]) . In each region, the deformations result in a vector RH problem with a 
sectionally analytic solution m^ a (z; to, t) which becomes unique given a generically true technical assumption 
discussed in Section [5j specifically Lemma 543 Here a stands for characters used to denote the asymptotic 
region (n, t) lies in ( e.g. a = cs for the collisionless shock region.) We often suppress the (n, independence of 
these vector-valued functions m^ a (z‘,n,t) = m^ a (z). 

We use the notation m Ki a(z) for unknown vector functions obtained through the deformation procedure 
where the integer k indicates how many deformations have been performed with k = 1 being the first defor¬ 
mation of RH Problem [3j The final deformation replaces the number k with the symbol ft. The subscript 
characters a. are used to denote the region as described above. For example, we will have the following 
sequence of deformations in the dispersive region: 

fh(z) i—mi )d (z) i—> m 2 ,d(z) i— m^z). 

The functions m K}a (z) are always related to m(z) via explicit transformations but they may not satisfy an 
RH problem with continous boundary values. Then m^ a (z) will always solve an RH problem with continuous 
boundary values and it is computed numerically. 

The jump matrix J(z;n,t ) in RH Problem [3] has terms that are highly oscillatory for most values of (to, t) 
and we need to control these oscillations in order to compute fh(z) (or, equivalently m(z)) accurately. To do 
so, we employ Deift-Zhou method of nonlinear steepest descent and examine the phase 0{z) that appears in 
these expressions. Solving 6'{z) = 0 for z , the stationary phase points of 9(z) are found to be 




= -™± 
t 


- -1 


Note that if n = t or n = —t, the stationary phase points coalesce at z = —1, or z = 1, respectively. 

Also, we present the results and deformations for the case n > 0 and t > 0. It is straightforward to obtain 
the solution for the case n < 0 and t > 0 by modifying the initial data and using n > 0 


Remark 4.1. If (a(t), b(t )) solves the Toda lattice with initial data a° and 6° and (a(t), b(t)) solves the Toda 
lattice with initial data a°_ n and — b°_ n+1 then a- n (t) = a n (t ) and b- n (t) = —6 n _i(i). An alternate approach 
would be to use the other reflection coefficient R-(z ) for n < 0 and let the stationary phase points lie in the 
right-half plane. 


We proceed with the details of the deformations used in each region. 


4.1. The Dispersive Region. In this region, the stationary phase points z^ 1 of the exponent 9{z) he on the 
unit circle, T. We set T, = {z: \z\ = 1, —1 < Rez < Rezo} and T = T \ S as shown in Figure [3] (this figure 
omits the contours D±). Note that since Re# = — Re#(2:), the curves Re#(z) = 0 are symmetric with 

respect to the mapping z ^ z~ x . 


Assume that we have performed the initial deformations discussed in Section 2.5 and obtained the vector¬ 
valued unknown fh(z',n,t). We then proceed with a deformation which will move the the oscillatory jumps 
along T into regions where the oscillatory terms decay exponentially. We define 

t ( z ) = 1 — -R(z)ii(2 -1 ) , 

and note that the jump matrix J(z m ,n,t) on T admits the following two factorizations: 

(24) J(z-,n,t) = M(z;n,t)P(z;n,t), 

where 


M(z; n, t ) = 


-R(i 


-1 


) e 

1 




and P(z; n, t ) = 


1 


R(z)d 


0(z;n,t) 
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Inu 



Figure 3. Sign of 8(z) and original jump contours for the RH Problem [I] in the dispersive 
region, n > 0. 


and 

(25) 

where 


L(z] n, t ) = 


1 

R(z)e e ( z ’ n ’ t '> 

K T 0) 


J(z;n,t) = L(z;n,t)D(z)U(z;n,t), 

/ -R(z~ v ) e - e ^ z ’ n ’ t '>\ 


, U(z;n,t ) = 


t(z) 

1 


and D(z) = 


t ( z ) 

0 


0 

l 

r(z) 


We use the MP-factorization (24) on T, and the LDC/-factorization ( |25| ) on £. M (for ‘minus’) will be 
deformed into the exterior (‘minus’ side) of the unit circle and P (for ‘plus’) will be deformed into the interior 
(‘plus’ side) of the unit circle. Then L is lower triangular and will be deformed into exterior of the unit circle, 
D is diagonal and will not be deformed, and U is upper triangular and will be deformed into interior of the unit 
circle. We employ these factorizations so that only one of or eappears in each matrix, which in turn 
makes it possible to obtain exponential decay in different regions of the complex plane. We introduce “ghost” 
contours, £±, deformed into ± side of £, and T± deformed into ± side of T. Note that these new contours 
pass locally along the directions of steepest descent for e± 9 ( z ’ n ’ t \ The first transformation in this region now 
follows. We define a new vector-valued function mi : d(z) based on the regions of the complex plane that 
emerge from this deformation, as shown in Figure |4](a). Note that mi i( j(z) satisfies the asymptotic symmetry 
condition and the quadratic normalization condition at infinity which are present in RH Problem [3j When 
j E K n ,t. we use the jumps Xj.± on Df, respectively, in order to turn exponential growth into exponential 
decay as t — > +oo. We use the jumps Yj- 1 on otherwise. 

More precisely, satisfies the following jump conditions: 

1 - 1 / 


m 1>d (z-,n,t)Q 
(z;n,t)Q~ 1 


m 


l,d 


(• z;n,t ) = < 


m 


l,d 

m i,d {z\n,t)Q- 1 

m Rd( Z ; n ’ t)Q~ l 
m~ A {z\n,t)Q- 1 

m i,d (z;n,t)Q~ 1 
m~ d (z-,n,t)Q 


-i 


[z)L(z; n, t)Q(z), 
z)D(z)Q(z), 
z)U(z-n,t)Q(z), 
[z)M(z; n, t)Q(z), 
z)P(z; n, t)Q(z), 
[z)X jj± (z-,n,t)Q(z), 
z)Y j) ±(z; n, t)Q(z), 


z E £_, 

z£S, 

2: E T_, 

z ^ r + , 

2 E Df, j E K nt , 


z E Df, j £ K, 


n.ti 


as 


seen in Figure Qb). Note that the definitions of T± and £± are given in the figure. 


Remark 4.2. The procedure by which the analyticity of an algebraic factorization of the jump matrix is 
exploited to modify the contours of an RH problem (cf. the transformation from m(z;n,t ) to m\^(z',n,t) 
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Figure 4. (a) Jump contours (blue) and matrices for RH Problem [3] with ‘ghost’ contours 
(dashed black), (b) Jump contours and matrices for This figure contains the definitions 

of r± and £±. 
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shown in Figure [4]) is referred to here as lensing. A full discussion of this can be found in [7J p. 191] although 
the term lensing does not appear there. 

Remark 4.3. Away from the points z^ 1 where the main contribution of the jump matrix is supported, 
the deformed contours E-t and r± are deformed to stay as far away as possible from T where J(z;n, t) is 
oscillatory for large t. More precisely, £± and r± pass from zg 11 locally in the directions of steepest descent 
of e =t0( 2 ;n,t). and away from z^ 1 , T,± and r± are chosen to be close to the inner and outer boundaries of 
the strip of analyticity of the reflection coefficient, while avoiding any intersection with the finitely many 
circles 1 < j < IV, or with the curves where R eQ(z;n,t) changes sign. This arrangement helps us gain 
sufficient exponential decay (to the identity matrix) in the jump matrices that are defined on the deformed 
contours. The analogous deformations in the Painleve, transition, collisionless shock, and soliton regions obey 
this principle. In the collisionless shock and transition regions, the points cr^ 1 and /5 ±:1 play the role of z^ 1 
for the arrangement of the deformed contours. 

Since generically i?(=tl) = —1 (see Appendix [E]) , the matrix D(z) has a singularity at 2 = — 1 (at z = 1 
in the case n < 0) and we need to remove this singularity. We also need the jump matrix to approach the 
identity to achieve accuracy for large values of the parameters. For these reasons we must remove the jump 
on the contour £ (see Figure [3]). This is achieved by solving a diagonal matrix RH problem with the jump 
contour £. We introduce the unique 2x2 matrix-valued function A(z) = A (z;n,t) that solves the diagonal 
RH problem: 


RH Problem 4. 

A + (z; n, t) = A _ (z; n, t)D(z), z£S, A(oo ;n,t) = I, 

(26) A (z; n, t) diag (|z + l| _1 ,|z + l|)=C>(l), z — 1, \z\ < 1, 

A (z; n, t) diag (|z + l|,|z + l| -1 )=0(l), z —> — 1, |z| > 1, 

such that A(z) is bounded for z in a neighborhood of zo, Zq 1 and the boundary values A ± (z), z G £, are not 
continuous only at zqjZq - 1 and — 1. 


It follows from classical theory that A(z) is a diagonal matrix with the property An(z) = 1/A22(z). The 
exact form of A(z), its properties, and a proof of the fact that it is unique can be found in the Appendix [A} 
Note that in general A(z) has singularities at the end points z = z^ 1 of the jump contour. To combat this 
issue we introduce circles around both Zq 1 , see Figure [5] We omit the conjugations by Q in Figure [H](b) and 
Figure [5](c) for the diagonal jump matrix D(z) since Q^(z), D(z), and Q(z) commute. 

We define m-z^z) as shown in Figure[5](b), where m2,d(z) = mi^(z) when no definition is specified; and we 
see that m2,d(~) satisfies the jump conditions that are presented in Figure [5jc). We apply the same procedure 
at Zq 1 . Along with the jump conditions in Figure [HJ m2,d(z) also satisfies an asymptotic symmetry condition 

(27) m 2l d(0) = m 2 ,d(oo) ^ Q( 0)A _1 (0) 

and the quadratic normalization condition at infinity present in RH Problem [3j Finally, we define m^d(z) 
by m^d^) = m 2 ,d( 2 )A _1 (z) and see that the vector-valued function mj ! d(z) satisfies the jump conditions 
shown graphically in Figure [6] This is the final deformation performed in the dispersive region. We have an 
important remark on boundary values of d(z) near z = — 1. 


Remark 4.4. Due to the singularity in L(z), D(z), and U(z) at z = —1 it is not immediately clear in what 
sense ra2,d(z) should satisfy the jump condition, or if a residue condition at z = — 1 is needed for 
Since we never solve for m 2 ,d(z) we can ignore this issue if we understand what conditions we need on m^d(z). 
From Proposition 2.2 and the behavior of A±(z) near z = — 1 we conclude that m(z)A _1 (z) (and hence 
m(z)A _1 (z)) is continuous up to T with jump 

(z)A)] 1 (z) = rh _ (z)Al 1 (z) [A_(z)L(z)Al 1 (z)A + (z)L r (z)Ai( 1 (z)] . 


m 


It then follows that both A^(z)L(z)Al 1 (z) and A+(z)t/(z)AT 1 (z) have analytic extensions to a strip that 
lies outside and inside the unit circle, respectively, and extend continuously up to the unit circle. Therefore, 
the jump contours and matrices of the vector problem (RH Problem [3]) can be deformed to those of m^d(z) 
leaving no singularity at z = — 1. 
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Figure 5. (a) ‘Ghost’ circles in preparation for the singularities of A ( 2 ), (b) Definitions of 
7712,(1 (*) near zq, (c) The jump contours and matrices for 7712 ,d( 2 ) near zq. 


One final detail to be covered is the radius of the circles we have placed near zq and z 0 1 . We follow the 
methodology put forth in [42] to determine this radius. Near zq we have 

6(z; n,t ) = 9(zo‘,n,t) - 2 ^ + (z - z 0 ) 2 + 0(z — zo) 3 - 

2 o 
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AQ~ 1 Y jt + QA ~ 1 




AQ-^j.-QA- 1 


X- 






AQ~ 1 Yj ) — QA ~ 1 


— 1 


_A C,/ ZA 





Figure 6. A zoomed view of the jump contours and matrices for the final RH problem in the 
dispersive region. 


We choose the radius of the circles to be proportional to r(n, t) := YF+^j so that f° r s := ( z — 2 o) / (cr(n, t )), 
| z — ^o| = |cr(n, t) | 

(28) 9(s ; n, t) := 9 (zq + scr(n, t); n, t)) = 6(zo',n, t ) — cs 2 (l + o(l)), 


where c is proportional to c 2 and accounts for the phase. Because 9(zo‘,n,t) is purely imaginary, e^ s;n ’^ is 
bounded if s is bounded. It also follows that for \n\/t < ci, in the dispersive region, we have r(n, t) = 0(t -1 / 2 ) 
so one may use r(n,t ) = cf -1 / 2 , in practice. 

The function 771^(2) satisfies a sectionally analytic RH problem with 


the jump conditions described in Figure [6j 

the asymptotic symmetry condition given in (27), and 

the quadratic normalization condition present in RH Problem [3} 


4.2. The Painleve Region. For? > 1 this region intersects with the soliton region defined below, and we 
use that deformation (see Section 4.5). For j < 1, the saddle points are coalescing at z = — 1 and this allows 
for a new deformation. Consider the arc £ that passes from z = — 1 and the two stationary phase points z J 1 
as shown in Figure [3J Set 2 = e l ( n ~ u ) and zo = eh 7r_a -’°) with cos(7r — 00 0) = — j. Thus z E £ if and only if 
|w| < |wo|, and wo —> 0 as t —> 00. Choose the branch cut (0, 00) for the logarithm. Then 9(z;n,t ) can be 
expressed in terms of w as: 


0(e*( 7r n, t) = 2 i (t sin(7r — w) + n(7r — w)) . 
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Note that > 1 — ct 2//3 implies 

| sin(vr - w)| < sinwo < V2ct~ 1/3 ^ 1 - ft” 2 / 3 = V2ct~ 1/3 (l - ^t~ 2/3 + ^t _4/3 + O (t -2 )) 
for large values of t. This together with t — |n| < ct 1//3 yields 

— 27rm = 2|tsin(7r — oj) — nw| < 2| sinwo|(i — |n|) < 2\/2c 2 + 0(t~ 2 / 3 ), 


for |w| < wo, which implies that the oscillations are controlled between the two stationary points. Therefore, 
the LZ)[/-factorization that was used in the dispersive region is not needed. Note that in this case the 
parametrix A( 2 ), which has unbounded behavior at 2 = —1 in the dispersive region, is not used. We perform 
a single deformation m(z) 1 —> m^p(z). Definition of the vector-valued function m^p(z) is given in Figure [7j[a). 
The jump contours and the jump matrices satisfied by m^p(z) are described in Figure[7j We note that m^p(z) 
takes continuous boundary values on its jump contour, and it is analytic on the arc of the jump contour where 
the jump matrix has been turned into the identity matrix, rn^p(z) satisfies a sectionally analytic RH problem 
with 

• the jump conditions described in Figure (7^b), 

• the asymptotic symmetry condition present in RH Problem [3] and 

• the quadratic normalization condition present in RH Problem [3} 


4.3. The Collisionless Shock Region. Recall that we use a deformation which involved A( 2 ) in the dis¬ 
persive region. The singularity at z = — 1 in the matrix D(z ) destroys the boundedness of the parametrbj^] 
A (z\n,t). As z —> —1, the matrices A(z)Q~ 1 (z)M(z)Q(z)A~ 1 (z) and A(z)Q~ 1 (z)P(z)Q(z)A^ 1 (z) are un¬ 
bounded and we cannot bridge the dispersive region and the Painleve region. By adjusting the constants 
that determine the asymptotic regions we can make the dispersive and Painleve regions overlap up to some 
finite t, but we wish to obtain a method which is stable for large values of t. To achieve this stability, we 
need to introduce additional deformations. The analogous region for the KdV equation has been introduced 
in [32] and the deformations were derived in [Qj. The asymptotic analysis of the solutions, the scaling, and 
the needed deformations for the Toda lattice in this region, to the best of our knowledge, are not present in 
the literature. 

As n increases in the dispersive region, the stationary phase points of e 8 ^ approach the singularity (2 = — 1 ) 
of the parametrix A(z). To prevent this, we replace the exponent 9(z;n,t ) by a so-called ^-function as was 
done for KdV in [13] (see also [32]). In what follows, we define the < 7 -function g{z) as the solution of an RH 
problem with the properties that mollify the unboundedness of D(z). Having done that, we introduce the 
needed deformations in this region. We leave the implementation details for solution of the RH problem given 


in (30) to Appendix C.3 In Appendix C.l we explicitly construct the ^-function, but it is more convenient 


to compute it numerically from the RH problem formulation given in (30). 

For n/t < 1, we let Ao G [—1, 0] and po G [0,1] denote the real and imaginary parts of the stationary phase 
point zq, respectively. Explicitly, 


(29) 


Ao 


n 

t 


and po 



For z 1 , Z 2 G T, z\ 7 ^ 22 , with 0 < arg Zj < 2n, j = 1, 2, define 

{z \, 22 ) arc = |e* e : min (arg z\, arg 2 2 } <0 < max{arg z\, arg z 2 } j 

oriented from z\ to 2 2 . Then [zi, 2 2 ] arc is defined to be the closure of ( 2 i, 2 2 ) arc . For a,j3 G T with —1 < 
Rea < Ao < Re/3 < 1, we label = [/3,a] arc , E c = [a, a -1 ] arc , and E/ = [a -1 , /3 -1 ] arc , as shown in 
Figure [8] 

Before describing the sequence of deformations used in this region, we proceed with the RH problem that 
determines the ^-function. Following the approach in [42] , we determine a and (3 on the unit circle so that 


o 

We use the term parametrix in a different way than is typical in the asymptotic analysis of RH problems. We use the term 
for any function that solves, or regularizes, any portion of the RH problem. 
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Figure 7. (a) Definition of m,jp(z) and the ‘ghost’ contours in preparation for the deforma¬ 
tion, (b) The jump contours and matrices for the final RH problem in the Painleve region with 
0 < n < t. 
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Figure 8. The cuts S u , £ c , and for 0 < n < t. 


there exists a function g(z) that satisfies the following properties, for some complex constants <5i and 82 


9 + (z) + 9 (z) = 


5i/t if 2 € T, u , 
—5i/t if z G £ i, 


(30) 


9 + {z) - g (z) = 6 2 /t, z <E £ c , 

g(z) - ^0(z) analytic in z for z 0 [/3,/3 _1 ] arc = S„US C U £/, 
g(z) is bounded at z = a ±l and z = P ±l , 
g(z) = \z — Ao log z + O (z _1 ) as z —> 00 , 

The constants S\ and 82 depend on a and (3 and, as we will see, they have the desired properties to eliminate 
the singularities. We leave the details of our method to solve this RH problem to Appendix C.3 Once g(z) 
is obtained, define the scalar function 


and construct the matrix function 
(31) 


0 ( 2 ) : = * (si*) - > 

^ Z) = { 0 e -*«) ’ 


which has the asymptotic behavior 

(j>(z) —> I as z —> 00 . 

Note that, for z £ T the jump condition satisfied by fh(z)(f>(z) is </>I 1 (z)Q _1 (z) J(z)Q(z)<f>+(z). In order to 
determine 8 \ and 82 we proceed as if Q(z) = I. It will be clear that this is sufficient. In this case, the jump 
condition (satisfied by the vector function m(z)</>(z)) on T is given explicitly by 


<t>- (z)J(z)</>+{z) = 


([1-R(z)R( Z - 1 )] (2) -R (z' 1 ) e -0(*)-0 + (*)-0 (O' 


^( z ) e 0 O)+s + (^)+s 0) 


=>-0 + O)+0 (z) 


because g(z) satisfies 


0 + (^)-0 (z) = t(g + (z)-g (z)) = 0, for z £ [/3, P ^ 
2g(z) = 2 tg(z) — 9(z) —> 0 as z —> 00 . 
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We write 


(32) <j>_ 1 (z)J(z)<f> + (z ) = 


'l- R(z)R(z~ 1 ) -R{z- 1 )e~ 2t ^ S 

R{z)e 2tfl ( 2 ) 1 

'[1 - R{z)R (z- 1 )] e <9 + (*)-9-(z)) _ R ( 2 -i) e -«i 


R{z)e‘ 


Si 


j(-9 + (z)+9 (z)) 


' [1 - R(z)R (z- 1 )] e 5 * -R (z- 1 ) e -*(ff + W+ff-W)' 
ii(z)e t ( 9+(2)+9 ' (0) ) e"* 2 

'[1 - R(z)R ( 2 - 1 )] e^W-s'W) -i? (^- 1 ) e 51 


R(z)e 


-Si 


j(-g+(z)+g (z)) 


'l - R(z)i? (z- 1 ) -i2 (z- 1 ) e“ w ®W' 
R{z)e 2t 1 


2 G (1, /3) arc , 

* G (^. «)arc > 

2 G ( a >« _1 )arc 

* G (a _1 ,/3 _1 ), 

(^ _1 , ^arc- 


Here 5i/t = g + (z) + g _ (z) for 2 G £ u , and = fl ,+ (^) — 9 ~{z) for z G S c . As can be seen in (32), this 
conjugation removes 0 (z; n, t) from the problem. 

We now present the initial deformation fh(z) 1 —> mi, cs (z) in the collisionless shock region. As in the 
dispersive region, we use the J(z) = L(z)D{z)U(z) factorization on H c and define mi jC8 (z), see Figure [9j(a). 
Here we used the lensing process (see Remark 4.2) to deform the RH problem. The jumps and contours near 
a -1 and /3 _1 are given in Figure tel b). What happens near a and (3 is clear by symmetry. 

We now perform our second deformation, mi iCS (,z) 1 —> m 2 ,cs(z), in this region. Define m 2 , cs (z) inside the 
circles centered at a and oT 1 as shown in Figure 10’a) and leave m 2 , cs (z) = rn\, cs (z) everywhere else. The 
jump conditions satisfied by m 2 ,cs{z) near the points a^ 1 and /3 -1 are shown in Figure |l0|(b). This deformation 
turned the jumps inside the circles surrounding a and a -1 to the identity jump, i.e. no jump. We remove 
0 (z; n, t ) from the problem using (j)(z), as discussed in the beginning of this section, by defining m^, cs {z) to be 


3±1 


m 3 , cs (z) = 


m 2 ,cs{z), inside the circles centered at cW 1 and /3 d 
m 2 ,cs(z)(j}(z), outside the circles centered at a^ 1 and /3 ±:L . 


The jump condition satisfied by m 3 , cs (z ) near the points a -1 and /3 _1 is shown in Figure 11 As in the 
dispersive region, the diagonal matrix D(z) has a singularity at z = —1 G E c since i?(±l) = —1 and this 
singularity has to be removed. 

We proceed with analyzing the jump matrix on E c in the limit 2:0 —> — 1 to determine the constants <5i and 
5 2 introduced in (30) (or in (32)) so that the singularity disappears. On E c = (a, o; -1 ) , the jump matrix is 

given by 


(33) 


4>_ 1 (z)D(z)4>+(z) = 


_ f[l-R(z)R(z~ 1 )]e! 


s 2 


0 


0 ([1 - R(z)R(z- 1 )\ e 52 ) 1 

Using R(— 1) = —1 and the analyticity of R(z) in a neighborhood around z = —1, we observe that 
(34) 1 — R(z)R (z _1 ) = v{z + 1 ) 2 + O ((z + l) 4 ) near z = — 1, 


for some constant v G C. So far, we have left a and (3 mostly arbitrary. It follows that (see Appendix C.l) 
the prescribed asymptotic behavior in (30) for g(z) as z —> 00 requires Rea + Re/3 = 2Ao, leaving us with 
single degree of freedom. Now consider the affine transformation, k = K(z), defined by 

'z~X 0 \ 


(35) 


K(z) = i 


P 0 


z{k) = K l (k) = Aq — ipok . 


Note that this transformation fixes the stationary phase points: K (z^ 1 ) = =pl, and the image of the contour 
(/3,/J -1 ) under the mapping K flattens as z$ —> —1 (see Figure 12). 

To remove the singularity of D(z) at z = —1, we need to obtain a parametrix ip(z) by solving the following 
diagonal RH problem: 
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(b) 

Figure 9. The initial deformation of the RH problem in the collisionless shock region, (a) 
Definition of m\ iCS (z) with its jump contours and matrices, (b) The initial jump contours and 
matrices near a -1 , (5~ l . 
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( a ) 



(b) 

Figure 10. (a) Definition of rn^^siz) near a -1 , /3 _1 , (b) The jump contours and matrices for 
m 2 ,cs {z) near ck -1 , /3 _1 . 


RH Problem 5. 

f + {z) = ip~(z)<j)Z 1 {z)D(z)(/) + (z), z <E S c , V’(oo) = 

ip(z) diag (\z + 1| _1 , \z + 1|) = 0(1), z- 1 , \z\ < 1 , 
ip(z) diag (|z + 11, \z + 1| _1 ) = 0(1), z— 1, |z| > 1, 

such that V’(^) is bounded for z in a neighborhood of the endpoints a, aT 1 and the boundary values ifr 1 {z ) 
are not continuous at a, a -1 and —1. 


In new variables (35), the jump condition in RH Problem [5] reads 


H + (k) = H~(k)D(k), k £ K (S c ), 


(36) 
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Figure 11. The jump contours and matrices for m 3, cs (~) near a , j3 1 



Figure 12. Cuts mapped under the affine mapping z k = K(z). 


where K~ 1 (k) = z(k), H{k ) = ip (K~ 1 ( k)), and D(k) = 4>Z 1 {K ~ 1 (k)) D(K ~ 1 (k)) 4>+ (K~ 1 ( k)). Let k* = 
K(—l) so that z{k) + 1 = —ipo(k — k*). We choose 62 to enforce p \^ 2 = 1 so that the (1, l)-entry of the 
diagonal jump matrix D(k) satisfies 

(37) [l — R (z(k)) R (z(fc) -1 )] e s-2 = v(k — k *) 2 + O [p^(k — fc*) 4 ) , near z =—1, k = k*, 


00 . In Appendix |C.1| it is 


hence removing, up to second order, the dependence on po- This indicates that D(k) remains bounded as 
Pq > 0 away from k* which ensures the boundedness of H, and hence of ip as t 
shown that a and j3 can (and should) be chosen so that 

- log Po 


(38) 


t 


[ 4 V(p~ a )(p- « _1 ) (p -P)(jp- /3 _1 ) + dp. 

h3 p “ 


See Appendix |C . 1| for the definition of the square root in ( |38[ ). Loosely speaking, the collisionless shock region 
is defined to be the region in the (n, t)-plane where ( [38] ) is solvable for a and (3 and a + 1 is not too small. 
This reasoning gives the asymptotic condition n = t — C 2 f 1 / 3 (logt) 2 / 3 . See Appendix 


C.2 


for more detail. 


Once is obtained (see Appendix 0, we conjugate the problem by ij)(z) as was done with A(z) in 
Section 4.1 Define m^ cs (z) by 


m $,cs{z) = 


n^3,cs(z)ip 3 (z), 
m3, cs (z), 


outside the circles centered at a and a 1 , 
inside the circles centered at a and a _1 . 
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The final deformation for this region and the RH problem satisfied by m^ cs (z) is shown in Figure 13 Similar 
to the case addressed in Remark 4.4 the jump contours and matrices of the vector problem (RH Problem [3]) 
can be deformed to those of m^ cs (z) leaving no singularity at 2 = —1, despite the fact that D[z ) and ^{z) 
are singular at z = —1. 





-i 

X 


0 

x 


'©'0 * 0 >O 




ip<l>- 1 Q- 1 X j , + Q<l> i/>~ 


0 © 


1 

X 




^<t> 1 



Figure 13. A zoomed view of the jump contours and matrices of the final deformation of the 
RH problem in the collisionless shock region. Note that 4>(z) and ip(z) commute. 


Finally, the choice of the radii of the circles round a, (3 , a _1 and f3~ 1 must be specified. It is easily seen 
from (pTOj ) that g'(z) vanishes as a square root at each of these points and g{z ) = a + b(z — c) 3 / 2 for c = a, 
/3, a -1 or /3 _1 and a, b depend on the choice of c. Following the arguments in (28) we choose the radius of 
these circles to be proportional to f~ 2//3 , of course, under the constraint that the circles should not intersect 
one another. 

The function m^ cs (z) satisfies a sectionally analytic RH problem with 

• the jump conditions described in Figure [l3j 

• the asymptotic symmetry condition 

"0 1 " 

1 0, 


(39) 


mj| iCS (0) = m(t iCS (oo) 


Q(O)0(O)^“ ± (O), 


and 
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• the quadratic normalization condition present in RH Problem [3] 


4.4. Transition Region. Similar to the case for the KdV equation (see [12]), the deformations in the colli¬ 
sionless shock region extends the values of (n, t) for which there exists a well-behaved RH problem beyond the 
dispersive region. However, this is not asymptotically reliable as we approach the Painleve region: as |n| — t 
decreases, a and or 1 approach the singularity of the parametrix ip(z) (see RH Problem[5]) at 2 = —1. To avoid 
this issue, we collapse the lensing on £ c = (a, a _1 ) arc that was introduced in the collisionless shock region (see 
Figure [9]) . Thus the LD 17-factorization of the jump matrix J{z\n,t ) is not used in this region. In order to 
maintain numerical accuracy, we choose a to ensure that the oscillations are controlled on [/3, /3 _1 ] ar[ .. The first 
deformation m{z) 1 —> mi^(z) we perform in this region is similar to the first deformation 777 .( 2 ) 1 —> mi tCS (z) 
in the collisionless shock region, but without the lensing on (a,a _1 ) arc . Definition of mi^(z) and the jump 
conditions it satisfies are given in Figure [14] 

Our second deformation involves conjugation by 4>(z) as in the collisionless shock region. We define ni 2 ,t(z) 
by ni 2 ,t(z) = mi )t (z)(j)(z), where 4>(z) is defined as in ([oT]) but we modify the definition g(z) below. The jump 


contours and the jump matrices for ni 2 ,t(z) near a 1 and /3 1 are presented in Figure 15 We will now show 


that collapsing the lensing on E c and conjugating by 4>(z) results in a well-behaved RH problem when the 
values (n, t ) lie in this region. In the analysis that follows, we omit the factors that come from conjugation 
by Q(z) to simplify the notation. As in the collisionless shock region, doing this has no effect on the result. 
Let n = t — t 1 / 3 r(t), where r(t) satisfies 

v(t) 

lim --= 0, and lim ?’(t) = 00 . 

t-yoo (logt) 2 / 3 roo w 

Given a positive bounded function f{n,t) we choose a and /3 by enforcing (recall that Rea + Re/3 = 2 Aq) 


(40) 


f(n,t) 


= —1 


i [ 4 V(P ~ «) (P - a -1 ) (P -P)ip~ /3 _1 ) 

.7-1 P 2 


dp. 


In light of this is equivalent to the conditions 
(41) 


t{g + (z)+g ( z))=if{n,t ) for 2G(/3,a) arc , 
t{g + (z) +g~(z)) = -if(n,t) for zg(^ 1 ,/T 1 ), 


By adjusting /, (40) can be solvecj^Jfor a since the right hand side is a monotone function of Rea under the 
constraint Re a + Re /3 = 2 Re Ao- 
Define h(n, t ) by 


h(n, t ) 


t 


[ 4 V(p - «) (p - « _i ) ( p -p)(p- / 3 _i ) 

Ja P- 


dp, 


and we have the following properties for g(z): 


(42) 


+ J if{n,t)/t if 2 € £ u , 

g W + s W = 

g + {z) - g~(z) = h(n, t)/t, z G S c , 

g{z) - ^0( z ) analytic in 2 for 2 ^ [/3,^ _1 ] arc = T, u U E c U E/, 

g(z) is bounded at 2 = a ±l and 2 = (d ±l , 
g(z) = \z — Aq log 2 + 0 ( ;2_1 ) as 2 00 , 


®In practice we use f{n,t ) = 2. Other choices may result in more efficient computations. 
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( a ) 



Figure 14. The initial deformation of the RH problem in the collisionless shock region, (a) 
Definition of mpfx) with its jump contours and matrices, (b) The initial jump contours and 
matrices near a~ 1 , /3 _1 . 
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Figure 15. The jump contours and matrices for m 2 ,t( z ) near a , (5 1 


Now, after applying the conjugation by 4>(z), again assuming for simplicity that Q(z) = I, as in the 
collisionless shock region, the jump matrix on T in this region is of the form 

( 1 -R{z)R{z~ 1 ) -R(z- 1 )e~ 2t ^ z y 
, R{z)e 2t9 ^ 1 


(43) <j>- 1 (z)J(z)<j> + (z) = * 


R(z) e ifM t 

' [1 - R{z)R (^~ 1 )] e h ("■*) -R (z- 1 ) 
R(z)e< 9+ W +9 ~W) e 

' [1 - R{z)R (^~ 1 )] e t{9 + { Z )-9-{z)) _ 

R(z)e ~ if ( re ’*) e 

'l -R(z)R(z~ 1 ) -R(z- 1 )e~ 2t9 ^ 
R(z)e 2t 1 


1’ 

z ^ (1) /-Oarc ! 

R (z -1 ) 

,t(-g + (z)+g~(z)) J 

> 2 e (0, a) arc , 

e -t(g + (z)+g~(z))\ 
—h(n,t) J ’ 

z E (a, a _1 ) arc , 

R (z- 1 ) e^M\ 
l (-g + (z)+g~(z)) J ’ 


)• 

Ze ^~^ 1 ) a rc- 


Note that (41), along with the fact that 

t\g + ( z ) + g~(z)\ <\f(n,t)\ for z E (a, a_1 ) arc , 

implies that oscillations in the off-diagonal entries of the jump matrix are controlled on (/3, /3 -1 ) . To analyze 

the situation concerning the diagonal entries, we find that t|g + (z) — g~(z) \ < \h(n,t)\ for z E (/3, a) arc . Using 
the change of variables z[k) = K~ 1 (k) given in (35), one can see that there exists a constant C > 1 such that 

1 < f(n,t) h{n, t) < c 

C ~ tpl tpl 

in this region, where po = Im zq as before. Now, note that 

3/2 


tpl ~ tV 8 1 - 


t — 


fV 3 r(f)\ 


= tV8t 1 r(t) 3/ ' 2 = \/8r(t) 3//2 —> oo, as t —> oo, 


by the assumptions on r(t ). This implies that 

tPo 


0 as t —> oo, 


and (40) is solvable for sufficiently large t. Furthermore, 

h(n, t ) ~ Ctpl —> oo as t —> oo, 
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which implies that the (2, 2)-entries of the jump matrix given in (32) all tend to 0 as t —> oo in this region. 
We are now left with the analysis of the (1, l)-entries of the jump matrix. We examine 

[1 -R(z)R(z~ 1 )]e h ^ on (/3,£ _1 ) arc , 

using the change of variables z(k) = K~ l {k). Observe that 

[1 - R{z{k))R{z{k)~ 1 )] e ft(n - f) = v(z{k) + l) 2 (l + 0(z(k) + 1 ) 2 )e^ = up 2 (l + O {pi)) e h 

uniformly in k for A = K(a ) and B = K{/3) bounded as t —> oo and K is defined in (J35J). Thus we are led to 
examine the behavior of p 2 e /i(n,t) f or i ar g e values of t > 0. Note that h(n,t) = 0{r{t y^), pi = 0(t~ 2 / 3 r(t)), 
and that for any c > 0 there exists T such that r(f ) 3 / 2 < clog(i) for t > T. There exists C\, C 2 > 0 

p 2 e h (n,t) < Cit -2/3 r ( t yC 2 clo g t < ^-2/3^0^ as t ^ 

for c chosen sufficiently small. This implies that the (1, l)-entries of the jump matrix all tend to zero. Therefore 
the entries of the jump matrix remain bounded and this gives us an asymptotically well-behaved RH problem 
without any lensing on [/?, /3 -1 ] arc - 

We now proceed with the final deformation in this region. We define 

m 2 ,t {z)^)~ 1 {z), inside the circles centered at a^ 1 and f3 ±l , 
rri 2 t t(z), otherwise. 


™u{ z ) = 


The jump contours and the jump matrices for the final RH problem for m^t(z) is given in Figure 16 The 
scaling of the circles around a, f3, a ^ 1 and f3 is the same as in the collisionless shock region: oc t -2 / 3 . 

The function m^ t {z) satisfies a sectionally analytic RH problem with 

• the jump conditions described in Figure [T 6 [ 

• the asymptotic symmetry condition given in ( |39| ), and 

• the quadratic normalization condition present in RH Problem [3j 


Remark 4.5. In the collisionless shock region when I := — log 

tpo 


0 , a « f3 and the g -function is essentially 
zero. This is the degeneration of the collisionless shock region to the dispersive region and hence the two 
regions overlap. The transition from the collisionless shock region to the transition region can be seen as 
t —> oo, when 


r A \]{<1- A ){q - B ) {q + A) {q + B) 


t = —i 


and A « 0 (see Appendix 


dq « —i 


i f \J(q — A) (q — B) (q + A) {q + B) dq, 
J B 


(A 0 - ipoq ) 2 

for the definition of A and B). This occurs when n = t — C\t 1 ^ 3 (logt) 2 ^ 3 as 


C.2 


Ci | 0 and it can been seen that the jump matrices in the transition region are regularized for r(t) = e(log t) 2 ^ 3 
for e sufficiently small. Thus the collisionless shock region and the transition region overlap. In the transition 
region as t —> oo we have 


f(n,t) _ \/(« - - B ) (<1 + A ) (? + B) 


tpl 


dq 


\j{q - A )(q - B ) (q + A ) {q + B ) dq- 


Ja (Ao - ipoq ) 2 jo 

Thus for f (n,t)/{tpl) sufficiently small, this equation is solvable for A. Then note that for n = t — cst 1 ^ 3 , 
tp 3 ~ ( 2 C 3 ) 3 / 2 . Thus choosing C 3 sufficiently large, we see that the Painleve region and the transition region 
overlap and the ^-function degenerates to zero as A —> 0 (or a —> — 1). In this way, our deformations can be 
seen to bridge all regions. An animation showing the deformations is given in the supplementary material. 


4.5. Soliton Region. In this region, we have |n| > t. Note that t = 0, \n\ > 0 is within the region. Therefore, 
the stationary phase points are no longer on the unit circle. Instead, for n > 0: 


n 

zo = -- + 



1 G (-1,0), 


and hence Zq 1 G (— 00 , — 1 ). Let A u = {re lU} : uj € [0, 2i r), r G (1 — u, 1 + i')}, v > 0, be the strip where 
R{z) is analytic. Note that such a strip exists as a consequence of our exponential decay assumption on 
(a^ — 1/2, 6 3 ). We have only one deformation to perform. We use the M P-factorization J = MP within A u 
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00 


— 1 

—X - 


00 


<l>- 1 Q- 1 Y :i , + Q4> 

0 

-X- 


00 


1 

-X- 




00 





PQ0 


Figure 16. A zoomed view of the jump contours and matrices of the final deformation of the 
RH problem (for i n the transition region. 


and deform T into two contours, T + and T_. As described in Remark 4.3| if zq and z 0 1 lie inside A v , T + 


and T_ pass through zq and Zq 1 , respectively, locally in the directions of steepest descent of e ±e ^ z ’ n,t \ Away 
from zq and Zq 1 , T + and are concentric circles that stay close to the inner and outer boundaries of A v , 
respectively. Once z^ 1 leave the strip A u deformed contours r± truncate to the concentric circles placed close 


to the boundaries of A v . Choices of T± are depicted in Figure 17 
As in the dispersive region, we use the jump 


M(z; n, t) 


^1 —R (z” 1 ) e ~ e ^ n ^ 


on the contour T_, and 

on the contour r+. We define the vector-valued function m^ s (z) in this region as given in Figure [XTjhi). 
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There is one final detail left to be covered concerning the signature of the real part of the exponent 6(z] nt). 
Define Co £ (—1,0) by R e6(z;n,t) = 0, that is, by 


(44) 


n 

t 


-i 


Co ~ Cp 

2 log | Co I 


Note that for n > t > 0, we have z o < Co (see [2TJ ). In light of the discussion in Remark 4.3, we arrange 
T ± C A v in a way that they do not intersect the curve given in (|44|). This ensures that the exponents in 


P and M have negative real parts on their domains. Consequently, P and M tend to the identity matrix 
exponentially fast as t — > oo. 

The strip of analyticity, jump matrices, and jump contours for the RH problem satisfied by m^ s {z) in this 
region are presented in Figure ]l7|(b) in the absence of poles {Cj}- When poles are present, we make sure 
that r± do not intersect with the circles, Dj, around each pole, j = 1,2,..., iV. We use the jumps 
that include and Yj ^ (defined in (21) and (22)) on D^ in presence of the poles precisely as described in 
Section ED 

The function m^ s (z) satisfies a sectionally analytic RH problem with 

• the jump conditions described in Figure [T7](b) , 

• the asymptotic symmetry condition present in RH Problem [3] and 

• the quadratic normalization condition present in RH Problem [3j 

5. Numerical solution of the deformed vector RH Problem 

In this section we describe, in detail, the methodology used to numerically approximate the solution m^ a (z) 
of the deformed vector RH problem and then produce the associated approximation of the solution of the 
Toda lattice. 


5.1. The numerical solution of the associated matrix RH problem. The asymptotic condition (12) 


is not convenient for numerical methods because it is nonlinear. But, as pointed out in [8], we can convert a 
vector RH problem such as RH Problem [3] to a (2 x 2)-matrix RH problem with the same jump conditions 
and a standard (linear) condition at oo provided this matrix problem has a solution. Then the solution of 
the vector RH problem can be reconstructed from the solution of the matrix problem. This reconstruction is 
discussed in the following section. 

Consider an RH problem 

(45) T + (z) = $~(z)J(z), zeT, $(oo) = /, 

which has smooth solutions (see Section 2.7 in m for the requisite conditions on J). We use [J;T] to refer 
to this matrix RH problem with the identity matrix condition at infinity. 

Definition 5.1. If a vector RH problem has the same jump condition as (45) then [J;T] is called the associated 


matrix RH problem. For example, RH Problem[8]below is the associated matrix RH problem for RH Problem]?] 

Given an oriented, piecewise-smooth contour F we define the Cauchy integral 

„ „ N If u(s) , 

Cru(z) = - / - ds. 

2m J v s — z 


It is well known that the operators defined by 


± 


Cp u(z) = (Cru(z)) 

are bounded operators from L 2 (T) to itself. Moreover, these operators satisfy the identity 
(46) C+ - Cf = /. 

If we assume that the solution to a matrix RH problem is of the form <f> = / + Cru, we can substitute this 
into the jump condition 4> + = 4> _ J and use this identity to obtain 


(47) 


C[G; r] := u - Cfu • (J - I) = J - I. 


This is a singular integral equation (SIE) for u. This motivates the following definition. 

Definition 5.2. The matrix RH problem [J; T] is said to be well-posed if C[J ; T] is invertible with a bounded 
inverse on L 2 (T) and J — I e L 2 (Y). 
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(a) 



(b) 

Figure 17. (a) Jump contour (blue) and matrices for the initial RH problem with ‘ghost’ 
contours (dashed black) that guide the deformation,(b) The jump contours and matrices for 
the RH problem satisfied by m^ s (z) in the soliton region. This figure contains the definitions 
of the contours T±. 
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This singular integral equation is critical in both the numerical and asymptotic solution of RH problems. 
A reader looking for a more in-depth discussion of RH problems should look to [2] for an introduction and 
|3 [3 07] for a more advanced discussion. For numerical solution of RH problems, we refer the reader to 
[31103 02] 03]. A comprehensive discussion of the inverse scattering transform can be found in [301. 

We also point out that the singular integral equation formulation we use differs from that of m ■ Our 
formulation has the benefit that the operator u i—)• C^u ■ (J — I) can be applied exactly to a chosen basis 
where the operator considered in m given by u i-a C v (u( J — /)) does not have this property. This seems to 
give a mild increase in the convergence rate. See E3 Chapter 2] for a comparison of the theory for these two 
formulations. 

Consider the contour T = (J" =1 Tj where each Tj is either a line segment or a circular arc. Thus we 
restrict to considering contours where a sequence of Mobius transformations Mi,..., M n are known such that 
Mfc([— 1,1]) = Tfc. Let P m , = {cos(j7r/m) : j = 0,1,..., m} be the Chebyshev points and let T m (x) denote the 
m th Chebyshev polynomial of the first kind. The points (J ? Mj(P nj ) are called the collocation points. Note 
that if a is an intersection point of a subset {Tj,,..., T Jm } of the contours Ti,..., F ? then it will be included 
777 times in this union. We include it m times by using the notation a + Oe* 0 A where 6j k is the angle at which 
Tj k leaves/approaches a. Additionally, / = g is used if and only if f(a) = g(a ) for all a E (J ■ Mj(F n .) that is 


not a point of self-intersection and lim^o /(a + ee lBjk ) = lim^o <7(a + ce 1,6 ^) if a is a point of intersection. 

A function 'L is said to satisfy the RH problem (45) at the collocation points (J ? Mj(F nj ) if ’ll has continuous 
boundary values and A J. The framework of Olver j23J implemented in [0H] (see also m) is designed 
to return a vector Vj of function values at the mapped points Mj (P nj .) (with directions attached at intersection 
points), so that the function U : T — > C* XJ defined piecewise by 


(48) 

(49) 

satisfies 


7 Ij 

U(z)\ Fj =Y j n l T l (M~\z)), 

i =0 

U(M( P n .)) = Vj, 


• I + C F U is a bounded function in C \ T, and 

• I + C F U satisfies the RH problem (45) exactly at Mj( P nj ). 

Here the coefficients a* in the definition of U in (48) are determined by the condition (49). 
We describe the method in more detail. Similar to before, substituting 


(50) 


4> ss / + C r U 


into the RH problem and using (46) gives a linear equation for U: 


(51) 


U — C^U ■ (J - /) = J-I. 


A closed-form expression for the Cauchy transform of the basis Tj(M ■ 1 (fe)) [45 . Section 4] (see also m 
allows the discretization of this linear equation by evaluating the Cauchy transform of the basis at the points 
Mj(F nj ). However, a modified definition for the Cauchy transform is required at the self-intersection (or 
junction) points points To (which are included in the collocation points Mj(P nj .)), at which the Cauchy 
transform of this basis is unbounded. By assuming that the computed U is in the class of functions for 
which I + C F U is bounded, we can define the bounded contribution of the Cauchy transform of each basis 
element at the points To- It can be shown, with appropriate assumptions on J (see the product condition, m 
Definition 3.8.3]) that the numerically calculated U must be in this class of functions. Therefore I + C F U will 
be bounded and satisfies the RH problem at Tq, hence at all points in Mj( P n .). 


We use (50) to show that if U E L 1 (T) then 


lim z($(z) — I) = — —^ [ U(t)dt, 
z ^oc 2m J p 


by the Dominated Convergence Theorem provided 2 is bounded away from T. The integral on the right-hand 
side can be computed using Clenshaw-Curtis quadrature. This relationship is needed in what follows to 
reconstruct the solution to the Toda lattice from the solution of the RH problem. Observe the complimentary 
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fact: if 

$(z) = A 1 (z)A 2 {z), Ai{z) = I + Ai^z^ 1 + 0(z~ 2 ), z->oo, 

then 


(52) $(z;n,t) = I + (Ai :l + A 2 ,i)z l +0(z 2 ), z —> oo. 

Another aspect of the numerical solution of RH problems is contour truncation. Note that if J(z*) = I then 
the linear system (51) at this point becomes U(z*) = 0. From this one can rigorously justify the removal of 
contours from the RH problem on which || J — iW^nL 00 * s small at the cost of a small error. This is discussed 
in more detail in Effl, Chapter 2] and see [4TJ for a discussion of implementing this idea. 


Remark 5.3. From the results in [25 J it follows that spectral convergence (i.e. convergence that is faster than 
n~ k for any k where n is the number of collocation points) can be verified a posteriori for a well-posed RH 
problem by checking that the norm of the inverse of the discretization of C[J;T] grows at most algebraically 
with respect to the number of collocation points. In the computations for this paper, we noticed at most 
logarithmic growth of the condition number for this collocation matrix, with a maximum on the order of 10 3 . 


5.2. Construct the solution of the deformed vector RH problem. With a method in hand to compute 
the solution of a matrix RH problem, normalized to be the identity matrix at infinity, we show that in order 
to solve the vector RH problem one can first solve the associated matrix RH problem and then take an 
appropriate linear combination of the rows of the solution of the associate matrix RH problem. This is the 
generic situation but there are some technicalities so we take care in the following developments. 

Let J : r — > C 2x2 where the contour F satisfies T -1 := { z : z E T} = — T and the minus sign refers to a 
reversal of orientation. Assume the symmetry condition 



The RH problems we want to solve are of the following form (compare with RH Problem [2]) which is assumed 
to be uniquely solvable: 


RH Problem 6. For an oriented contour T, we seek a function v: C\T — >• C lx2 that is sectionally analytic, 
continuous up to T and satisfies: 

• the jump condition: 

(54) v + (z) = v~(z)J(z), zeT, 

• the symmetry condition: 

(55) v (z) = v(£ _1 ) , 

• and the normalization condition: 


(56) 


v(oo) = lim v(z) = [vi V2) , v\ ■ V 2 = 1 , v\ > 0 . 


In the background, throughout all of our deformations of RH Problem [2] is a problem of the form of RH 
Problem [6j But we do not preserve the symmetry condition through the deformations, mainly for convenience 
and ease of numerical implementation. We introduce the notion of a non-singular deformation to encapsulate 
this: 


Definition 5.4. A vector or matrix function v is a non-singular deformation of a sectionally analytic function 
v : C \ r —> C^ 2 , for j = 1 or j = 2, with continuous boundary values if there exists a sectionally analytic 
matrix function H : C \ T 7 —> C 2x2 , det H(z) = 1, also with continuous boundary values, such that 

(57) lim H(z) =diag(ci,c 2 ) 

z—t OO 

exists and 

v(z) = v(z)H(z), ze C\f, f:=(rur'). 
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So, if v : C \ f -A & x2 , j = 1 or j = 2 , is a non-singular deformation of v(z) then for 

(58) J(z) := Hl\z)J(z)H + (z ), 
v(z) satisfies: 

RH Problem 7. For an oriented contour T, we seek a function v: C\T —> C lx2 that is sectionally analytic, 
continuous up to T and satisfies: 

• the jump condition: 

(59) v + (z) = v~(z)J(z), z£ f, 

• the asymptotic symmetry condition: 

(60) v (0) H~ 1 { 0) = {5(oo)H _ 1 (oo) ^ ^ , 

• and the normalization condition: 

(61) u(oo)-ff - 1 (oo) = (ui V 2 ) , ui ■ V 2 = 1, v\ > 0 . 

A solution of RH Problem [ 6 ] clearly produces a solution of RH Problem [7] and so RH Problem [7] is solvable 
(we assume RH Problem plis always uniquely solvable). Since RH Problem [7| turns out to be a bit more 
numerically tractable, we want to know when the solution of RH Problem [7] is unique. So, consider the 
associated matrix RH problem: 


RH Problem 8. For a bounded, oriented contour T, bounded away from the origin, we seek a function 
'F: C \ r —> C 2x2 that is sectionally analytic, continuous up to T and satisfies: 

• the jump condition: 

V + (z) = V-(z)J(z), ze f, 

• and the normalization condition: 

’F(oc) = I. 

Lemma 5.5. Suppose that RH Problem [5| is uniquely solvable. Then RH Problem [7] is uniquely solvable if 
and only if RH Problem [6] is uniquely solvable. 

Proof. We begin with a straightforward calculation. Let v(z) be a solution of RH Problem [7| Define v(z) = 
v(z)H~ l (z), a solution of RH Problem [ 6 ] with an asymptotic symmetry condition (not the global symmetry 
condition (55)): 


(62) 

Then consider 


D(0) = D(oo) 


v(z) := v(z x ) 


0 1 
1 0 


0 1 
1 0 


We see from (53) 


v + {z) = v (z 1 ) 


= v + {z- l )J-\z~ 1 ) 


= v ( z)J(z ). 


Therefore v(z) satisfies the jump condition in RH Problem [ 6 ] with the asymptotic symmetry condition (62). 
Now, assume RH Problem [7] is uniquely solvable, with solution v(z). Define k(z) by 

"0 1 " 


k(z) = \v(z)H-'(z) + i itz-'JirV 1 ) 


1 0 


It follows that k(z) satisfies (56) and therefore k(z) is a solution of RH Problem [ 6 j It is clear that any solution 
of RH Problem [ 6 ] gives a solution, via // (z). of RH Problem [7] and so RH Problem [ 6 ] is uniquely solvable. 

Assume RH Problem [ 6 ] is uniquely solvable with solution v(z). Then v(z)H(z) is clearly a solution of RH 
Problem [7J Assume then that k(z) is a new solution to RH Problem [7] with k(z) 7 ^ v(z)H(z). Define 

r 


v(z) = \k{z)H-\z) + hiz-^H-^z- 1 ) 


1 0 






NUMERICAL INVERSE SCATTERING FOR THE TODA LATTICE 


35 


which is a solution of RH Problem [ 6 ] and by uniqueness v(z) = v(z). From (60) 

&( 0 ).ff _ 1 ( 0 ) = k( oo)Ff~ 1 (oo) 

so that 


0 1 
1 0 


v(oo) = -k(oo)H 1 (oo) + -k(0)H 1 (0) 


0 1 
1 0 


= k(oo)H (oo) = v(oo). 


So, let w(z) ■= (v(z) — k{z)H 1 (z))H(z) = v(z)H(z) — k(z), the difference of two solutions of RH Problem[7j 
It follows that w(z) satisfies the jump condition ( [59] ) with 

w(oo) = (0 0) . 

From the uniqueness of solutions of RH Problem [ 8 j w is identically zero. □ 

Given RH Problem [7] we call RH Problem [ 8 ] the associated matrix RH problem. The goal is to compute 
the coefficients in the expansion 

v(z) = (vi + Vl,lZ V 2 + V2,iz) + 0(z 2 ), z -> 0. 

Since we know that RH Problem [7] and RH Problem [ 6 ] are equivalent if RH Problem [ 8 ] is uniquely solvable, 
we use the following: 

Lemma 5.6 (From matrix solution to vector solution). Assume RH Problems [#| and[ 6 | have unique solutions 
T and v, respectively. Then the left nullspace of T := ’P(0)IL~ 1 (0) — H~ 1 ( oo) ^ is one dimensional and 

v(z) : = Y^(z)H~ 1 (z) 

solves RH Problem^ where Y is the left null vector o/T chosen so that X := YH~ 1 ( oo) satisfies X\ > 0 and 
Y-Y 2 = 1. 

Proof. Let v(z) be the unique solution of RH Problem [g] and we must show that v(z) = Y^>(z)H~ 1 (z) for 
a unique vector Y. As we assume RH Problem [ 8 ] is uniquely solvable with solution T, from Lemma |5.5| it 
suffices to enforce the asymptotic symmetry condition ( |60[ ) because the unique solutions to RH Problem [7] 
and RH Problem | 6 ] coincide. By enforcing ( |60| ) 

(63) YV{0)H~ 1 (0) = Y'lt(oo)H- 1 (oo)ft J 

Then Y, if it exists, must be in the left nullspace of \k := 'F(0)i7 _ 1 (0) — H~ l { oo) ^ ^. Now, to see that 

we can find a left null vector that can be normalized by X := YH~ 1 ( oo), X\ > 0 and X\ ■ X 2 = 1. Consider 
the matrix function 

._ ( 


: \v{z)H{z) / 

where subscript refers to the first row. Because neither v\ nor v 2 in v(oo) = (ui v 2 ) can vanish and 
'Fi(oo) = (l 0), I\( oo) is invertible by ( [57] ). We set Z(z) = K~ 1 (oo)K(z) so that Z(oo) = I and by 
uniqueness \k = Z. It suffices to take Y = (0 l) K(oo) and Y must exist. 

Now, assume there is another left null vector Y of T that is not a multiple of Y. Then T = 0 and Y\ ■ Y 2 >0 
so Y = (Yf 0) is linearly independent of Y. Let a > 0 and /3 > 0, and consider 

Ti7 _ 1 ((X)) = (3{Y + aY)H~ 1 ( oo) = (/3(Ti + aYf)/c cfdYf) , H~ 1 ( oo) = diag(l/c, c). 

As Yi • Y 2 = 1, Y\ = Y 2 * /1 T" 2 1 2 so that argld = arg Yf. Then, choosing a > 0 sufficiently small 

U >0 ^ W + a r 2 -) > o 
c c 

Then by choosing (3 so that Y\ -Y 2 = 1 we find that 

v(z) = Y'P(z)H~ 1 (z), 
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is another solution of RH Problem [ 6 ] with 'O(oo) ^ 'u(oo), violating uniqueness. Therefore Y is uniquely 
defined. □ 

There is a subtlety here that will become apparent as we proceed. We often know that the non-singular 
deformation carried through by H(z ) exists but we will not want to compute it. We will also encounter 
singular deformations, i.e. we will multiply the solution of RH Problem [2] by matrix functions that have 
singularities. Assume that 

v(z) = v{z)H{z) = fh(z)S(z)H(z), 

where S(z) is a, possibly singular, matrix function with det5(*) = 1. The product S(z )H( z ) is something 
we will know explicitly and be able to compute (i.e., S(z)H(z ) = A~ 1 (z)Q(z) in Section 4.1). Because m(z) 
and v(z) both satisfy the symmetry condition it follows that 


S(z) = 


0 1 
1 0 


S(z 


-i\ 


0 1 
1 0 


0 1 
1 0 


From (63), assuming S(z) is analytic in a neighborhood of z = 0, 

i/(0) = TT(0)^ 1 (0) = Y'Y(oo)H~ 1 (oo' / 

= YV(oo)H~ 1 ( 00 ) 

m(0) = YVfflH-^S-^O) = ^(ooJirHoo)^ 1 ^) 


1 ils-'o). 


0 1 
1 0 


To compute Y it suffices to know the product S(z)H(z) at infinity and at zero. Define P(z) = S(z)H(z ), 
T 1 := lim 2 _ s . 00 z(^>(z) — I) and P 1 := lirn^oo z(P(z) — P( 00 )). Then 


m(z) = Y^(z)P~ 1 (z) = YV(z- l )p-\z- L ) 


-In 


0 1 
1 0 


= [y^(oo)p- 1 (oo) + [zY T 1 - p-^oo^p-^oo)]] ^ ^ + 0(z 2 ), z -> 0 . 

One should think of RH Problem [ 6 ] as being an abstraction of RH Problem [2j Our basic assumption is that 
the associated matrix RH problem to a non-singular deformation of a problem of the form of RH Problem [ 6 ] is 
uniquely solvable. This assumption is not violated in practice but one cannot rule out exceptional cases, see 
Remark 6.1 To compute the solution of a vector RH problem with the normalization (61), we first deform 
the vector problem at hand, then solve the associated 2x2 matrix RH problem and then use that matrix 
solution to construct the solution of the vector RH problem. This process fails only if the solution of the 
associated matrix RH problem fails to exist. 


5.3. Extracting the solution of the Toda lattice. The numerical procedure here returns an approximation 
of m^ a (z', n, t ) = Y^f(z) for a matrix-valued function ^( 2 ) and a row vector Y. In a neighborhood of infinity 
m(z;n,t) = m$ j 0 l (z;n,t)Z(z) = Y''S>(z)Z(z) for a (locally) analytic function Z(z) such that Z( 00 ) = /. From 
( [52] ) with 

'P(z) = I + ^>iz~ 1 + 0(z^ 1 ), Z(z) = I + Z\z _1 + O(z ^), z -A 00 , 

we have 


m(z;n,t) = Y(I + (Ti + Zi)z ) + 0(z ). 


This combined with Lemma 2.3 and (23) is enough to compute the solution of the Toda lattice. 

Remark 5.7. In choosing Z( 00 ) = I, which simplifies this calculation, we cannot perserve the symmetry 


condition m(z; n, t ) = m(z 1 ]n,t ) 


0 1 
1 0 
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6. Numerical results 


6.1. Direct scattering. In this section we present some numerical results on the computation of the scat¬ 
tering data. We study two choices of initial data in detail: 

(TS) A choice of initial data giving rise to two solitons (TS) is 


a n (0) = ^ + ^ne~ n \ 

K (0) = ^ sech(n). 

(NS) A choice of non-solitonic (NS) ( i.e ., a pp (L) = 0) initial data is 

. . 1 1 _2 

«n(0) = - - , 

bn(0 ) = sech(n). 


The reflection coefficient on T is shown in Figure [l 8(a) for TS initial data and in Figure 18(c)| for NS initial 
data. With K in Section 3.1 sufficiently large (K = 30 is sufficient), accuracy is guaranteed. In the case of 
the TS data, we find 


(Ci.Ca) ~ (0.596142,-0.704859), 
( 71 , 72 ) ~ (3.25791,1.43054). 



Figure 18. Numerical computation of the reflection coefficient R(z ) (Rei2(z) solid curve, 
Imfi(z) dashed graph) with initial data where (a) two eigenvalues are present (TS) and (c) 
fjpp(L) is empty (NS). We plot 1 — |-R(z)| 2 for TS and NS initial data in (b) and (d), respectively. 
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6.2. Inverse scattering. In this section we present numerical results for the computation of the inverse 
scattering transform. These results are of three flavors: 

• Example solution plots, 

• error analysis, and 

• numerical asymptotics. 


6.2.1. Example solutions. Here we present plots of the solution of the Toda lattice with both TS and NS 
initial data. See Figures 19 and 20 for plots of a n (t) and b n {t) when t = 100,1000 in the case of TS initial 
data. Two solitons traveling in opposite directions are clearly visible. Indeed, this is anticipated because 
and C 2 have opposite signs. See Figures 21 and 22 for plots of a n (t) and b n (t) when t = 100,1000,10000 in 
the case of NS initial data. As stated above, no solitons are present in the solution and the high oscillation 
in the solution is apparent, especially at t = 10000. 



e 

e 


-1000 -500 0 500 1000 

n 

(b) 



Figure 19. Solution a n (t ) obtained by numerical inverse scattering (a) at t = 100, (b) at 
t = 1000. The horizontal axis denotes the spatial parameter n E Z. The initial data produces 
two eigenvalues on opposite sides of the a.c.-spectrum giving rise to two solitons traveling in 
opposite directions. 

6.2.2. Error analysis. To examine the accuracy of our numerical inverse scattering transform (1ST) a posteriori 
we compare it with a naive time-stepping method. A more detailed description of a related method can be 
found in [5j. Here we just use out-of-the-box Runge-Kutta 4. Fix K > 0 and consider the Toda lattice with 
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Figure 20. Solution b n {t ) obtained by numerical inverse scattering (a) at t = 100, (b) at 
t = 1000. The horizontal axis denotes the spatial parameter n E Z. The initial data produces 
two eigenvalues on opposite sides of the a.c.-spectrum giving rise to two solitons traveling in 
opposite directions. 

Dirichlet boundary conditions: a-x(t ) = ax(t) = 1/2, b_x = bx = 0. Here K is chosen sufficiently large 
so that the solution remains flat at the boundary ±K for all times simulated. If t < T, it suffices to take 
K = cT for some c > 1 which is larger than the speed of the fastest soliton present in the solution. This is a 
finite-dimensional system of ODEs and can be integrated in time using the fourth-order Runge-Kutta method. 
We compare the time-stepped solution at t = 30 with A t = 10“ 5 with the solution computed via the inverse 
scattering transform in Figure [24] As the number of collocation points in increased, the numerical inverse 
scattering solution converges exponentially to the true solution, with these errors saturating at approximately 
10“ 11 . It is reasonable to expect that at this point the numerical 1ST gives a more accurate solution than the 
time-stepping method. Furthermore, on a standard laptop it takes ~ 6 seconds to compute the solution at 
t = 30, n = 30 with 720 collocation points using the 1ST and ~ 2x 10 5 seconds with the naive time-stepping 
methocp^l implemented in Mathematica. Presumably, by using more efficient integrators and software packages 
this time can be reduced by at least an order of magnitude but when computing at sufficiently long times, the 
numerical 1ST is guaranteed to have a shorter runtime. This comparison is pessimistic as contour truncation 

should be noted that the time-stepping method produces an approximation of the entire solution profile in this time while 
the numerical 1ST gives the solution at only one point. Even so, it would take « 1.2 x 10 4 second to compute the entire solution 
profile at this rate with the numerical 1ST. 


n 


(a) 



n 

(b) 
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n 
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(c) 


Figure 21. Solution a n (t) obtained by numerical inverse scattering (a) at t = 100, (b) at 
t = 1000, (b) at t = 10000. The horizontal axis denotes the spatial parameter n € Z. The 
initial data produces no eigenvalues and therefore no solitons are present in the solution. 
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(c) 


Figure 22. Solution b n (t) obtained by numerical inverse scattering (a) at t = 100, (b) at 
t = 1000, (b) at t = 10000. The horizontal axis denotes the spatial parameter n 6 Z. The 
initial data produces no eigenvalues and therefore no solitons are present in the solution. 









42 


D. BILMAN AND T. TROGDON 


(see Figure 25 and the next paragraph) reduces 1ST computation times as t increases while time-stepping 
methods see their complexity increase. 

We emphasize that the number of collocation points required to solve the Toda lattice using the numerical 
inverse scattering transform to a given accuracy is typically decreasing with respect to t. This is because 
the deformations performed on the original RH problem force the jump matrices on some contours to tend 
exponentially fast to the identity matrix as t increases. Thus, after truncation, discussed in Section [6j fewer 
contours need to be discretized in the RH problem resulting in fewer collocation points. To see this in action, 
using the contour truncation algorithm described in [4lj . see Figure [25j Thus the errors seen in Figure [24] are 
pessimistic for large values of t. 


Remark 6.1. Given a uniquely solvable vector RH problem such as RH Problem [7j it does not follow that the 
associated matrix solution exists. It is well-known (see [45] . for example) that the operator C[J;T] associated 
with RH Problem [ 2 ] is Fredholm with index zero. Furthermore C[J -1 ; r]C[J; T] = I + I\ where K is a compact 
operator. For fixed n, K = K(t) is analytic in t and then by the analytic Fredholm theorem C[J;T] is either 
never invertible or invertible on the compliment of a discrete set values of t. From the work of [8], it can be 
deduced that for each fixed n, there exists t*(n) > 0 such that if t > t*(n ) then C [,/; T] is invertible. Therefore, 
we know that there is only a discrete set of possible values of (n, t ) where C[J\ T] may fail to have an inverse. 
This is discussed in Appendix |D| 

To investigate the possibility of encountering a point (n, t ) we perform the following computation. Fix t > 0 
and vary n. For each value of (n, t) we plot the smallest singular value of the discretization of C[J;T] found 
using the numerical method described in this section. It is known that C[J;T] is always invertible when no 
solitons are present — when no poles are present in RH Problem [l] In Figure [23] we perform this experiment 
for initial data with and without solitons. Because n is discrete, and we will only ever evaluate at discrete 
times, we do not ever expect to encounter a singular operator C[J;T]. 


Two Solitons No Solitons 



(a) (b) 

Figure 23. Plotting singular values for t = 0. (a) A plot of the singular values of the 
discretization of the operator C[J; r]from RH Problem [ 2 ] when the initial data produces two 
solitons. (b) A plot of the singular values of the discretization of the operator C[J; T] from RH 
Problem [2] when the initial data produces no solitons. In comparing these two plots, it is clear 
that the smallest singular value abruptly approaches zero, but due to the discrete nature, a 
zero singular value is not encountered. In panel (a), the minimum at n = 1 is ~ 0.0013. For 
the exact form of the initial data see (NS) and (TS) initial data in Section [6j 


6.2.3. Numerical asymptotics. We have seen that the accuracy of the numerical method is easily verified for 
short/moderate times by the comparison with time stepping methods. Thus in terms of accuracy, the method 
will also out-perform the asymptotic formulae ( [341 [33l fl9 j) for long-time. Thus, with a numerical inverse 
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Figure 24. The error between a solution of the Toda lattice computed using time stepping 
and numerical inverse scattering. The error shown is the maximum of the error computed 
at the points (n,t) = (—51, 30), (—31, 30), (—11,30) plotted versus the number of collocation 
points used in the numerical solution of the associated RH problem. 



scattering transform the numerical evaluation of asymptotic formulae (which is truly a non-trivial task) is no 
longer necessary in many cases. To demonstrate this, we show the long-time behavior of the solution with NS 
initial data in the dispersive, Painleve and collisionless shock regions. We call such computations numerical 
asymptotics. 


Dispersive region. 

To show the solution in the dispersive region we let n depend on t through n = \7t/ 10J where |_ - J 
represents the integer part. See Figure [26] for a plot of the solution into the dispersive region. 
Painleve region. 

To show the solution in the Painleve region we let n depend on t through n = \t — t 1//3 J. See Figure [27] 
for a plot of the solution into the Painleve region. 

Collisionless shock region. 

To show the solution in the collisionless shock region we let n depend on t through n=\t — 
t 1//3 (log t) 2 / 3 J. See Figure 28 for a plot of the solution into the collisionless shock region. 


Appendix A. Solving the singular diagonal RH problems 


There are two diagonal RH problems that must be solved, and computed numerically, in our deformation 
procedures. The first is A(z) in RH Problem [4] and the second is 'ip(z) in RH Problem [5j Note that because 
g + (z ) — g~(z ) is constant for z G (a, a _1 ) arc , these problems are both of the form: 

RH Problem 9. 

X+ ( z ) = X~{z) 2 e (7-T' 1 )^, M = 1, Im7 > 0, 

X(po) = /, 

X(z) diag (|z + l| -1 ,|z + l|)=C ) (l), z -A — 1, |z| < 1, 

X{z) diag (|z + 11, |z + lj^ 1 ) = 0( 1), 2 -A —1, |z| > 1, 

for a positive real constant c such that X(z) is bounded for z in a neighborhood of 7, and the boundary 
values T =t (2) are not continuous at and —1. 


It follows from classical theory that X \2 = X 21 = 0, X\ \ (z) = 1 jXyiiz). From this 


Xn(z) = exp 


27TZ J( 7;7 -l) a 


s — z I 


Furthermore, because log(cr(s)) is real valued it follows that X\ 1 (z) is bounded near 7 and 7 1 , see [23]. 
The difficulty here is evaluating X\\{z) numerically because logr(s) is singular at s = -1 S (7)7 -1 ) • We 
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Figure 25. Contours used in numerical inverse scattering for the Toda lattice in the dispersive 
region with TS initial data. This figure shows how contour truncation localizes the contours 
as t -» oo. (a) The contours for (n,t) = (—11,30) with 28 total contours, (b) The contours 
for (n,t) = (—110,300) with 24 total contours (soliton contours have been truncated), (c) The 
contours for (n,t) = ( — 1100,3000) with 14 total contours. 


regularize the integrand by considering^] 


t(z) = 


^ 0 ) 

(z — z~ x ) 2 


>0, ZG (7,7 ') arc , 


the amplitude of the initial data increases, f(±l) approaches zero. Round-off error that is amplified by forming the ratio 
may become significant and degrade the accuracy obtained in evaluating f. In this case, once should look for a different method 
to compute X\\. 
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Figure 26. Numerical asymptotics in the dispersive region: n = |_Tt/10j. Such computations 
are accurate for arbitrarily large t. 




101000 


105500 

t 


110000 


Figure 27. Numerical asymptotics in the Painleve region: n = [t — f 1 / 3 ] . Such computations 
are accurate for arbitrarily large t 
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Figure 28. Numerical asymptotics in the collisionless shock region: n = [t — t 1 / 3 (logi) 2 // 3 J. 
Such computations are accurate for arbitrarily large t. 


so that the Cauchy integral of the smooth function logf(z) along ( 7,7 1 ) arc can be computed with the 
methods referenced in Section [ 6 j Consider 

1 r log(s + l) ^ | 1 r log(l-s _1 ) 

s ~ z tt* 7( 7 , 7 - 1 ) a 


x 0;7) : = [ 

2vr* 2(7,7 Pa 


log(-(s-s 1 ) 2 ) d? = J_ 


1 

2ni 


(7-i) a 


s — z 
in 1 

~ dS+ 2ii 


-ds 


s — z 


ITT 


ds. 


h-hT-Pa 


s — z 


To verify the second equality, we choose the branch cut of log(l — z _1 ) to be on [0,1] with it being real-valued 
for z < 0. We then choose the branch cut of log(l + z) to be on (— 00 ,—1] with it being real valued for 
z > —1. Note that these are just the composition of the principal logarithm with z 1 —> 1 — z^ 1 and z 1 —y 1 + z, 
respectively. With these definitions, we must show that 


21og(z + 1) + 21og(l — z x ) — 7 risignlm£ = log (—(z — z 1 ) 2 ), \z\ = 1. 
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For Im z > 0, log (2 +1) + log(l — z _1 ) = log (z — z^ 1 ) and z — z^ 1 is purely imaginary with Im (z — z -1 ) > 0 so 
that 2 log (2 + 1) + 2 log(l — z~ l ) — ni signlmz = log(— (z — 2 V 1 ) 2 ). A similar calculation follows for Im z < 0. 

Furthermore, signlmz is piecewise-smooth on ( 7 , 7 _1 ) and log(l — z~ l ) is smooth so that their Cauchy 
integrals can be, again, computed with the methods referenced in Section [ 6 j Thus it remains to calculate by 
Cauchy’s integral formula 


(64) 


1 f 2 log(s + 1) ^ _ 1 /' 2 log(s + 1) ^ I 0, if 2 G T+, 

2m ./( 7 ) 7 -i) arc s-z 2m J^ s- z [21og(2 + 1), if z e D y , 


where J 7 is the vertical line connecting 7 and 7 1 with downward orientation and D 7 is the region enclosed 
by I 7 and ( 7 , 7 -1 ) ar • It turns out that 


2ni 


s — z 


2m V 


- log(l + 7 ) log 


z — 7 
1 + 2 


+Li2 


'1 + 7 


-1 ■ 


1 + ^ 


— Li2 


A + 7 n 
1 + z, 


where Li 2 is the dilogarithm function, see |2~4[ Section 25.12], with appropriately chosen branch cuts. Stock 
special function routines in Mathematica allow this function to be computed accurately. Then 


Pd 1 1 ( 2 ) = exp 


, ^ J ( 7 , 7 1 ) a 


log(cf(s)) 


s — z 


ds + l(z; 7 ) 


Finally, we examine the singularities of +’ 11 ( 2 ;). First, we note that if f(s) is a continuously differentiable, 
real-valued function on a contour C oriented from a to b then 


_L f J^l ds= -L/(s)log(s - z) 

2m Jq s — z 2m 


~ 7^7 f f\s) log(s - z)ds. 
J c 


As 2 -A b the only unbounded term is (27 xi) 1 f(b) log (6 — 2 ) and when exponentiated, 

exp [(2m)~ l f(b) log (b - 2 )] = (b - z )D b )/ 2 ™, 


is bounded as 2 — > b because f(b) is real. Therefore + 11 ( 2 ) is bounded near zq, Zq . We now must investigate 
the singularity at 2 = —1. The only singularity can come from I(z; 7 ). When exponentiated, (64) contributes 
a second-order zero as 2 —> — 1 from inside the unit circle. When exponentiated, the quantity 


1 

2ni 




1 f in 

2ni ]<_ 1 7 -i) s-z 

V ’ 1 /arc 


produces a simple pole at 2 = —1. Hence, we find that + 11 ( 2 ) has a simple pole as 2 -A —1 from outside 
the unit circle and a simple zero as 2 A —1 from inside the unit circle. It is then clear that X(z) does 
indeed satisfy the conditions set forth in RH Problem[9| Now assume y(z) is another solution. It follows that 
y(z)X~ 1 (z) has (possibly) isolated singularities at 2 = —1, 20 , Zq 1 but, the singularities must be bounded and 
hence this product is entire. By the asymptotic condition, y(z) = X(z) and X(z) is the unique solution. 


Appendix B. On computing eigenvalues of L 

In this section we present a case where we fail to capture all of the eigenvalues of L numerically by using 
conventional eigenvalue algorithms on finite, K x K truncations Lk of the doubly-infinite Jacobi matrix L. 
For 

KiL) < A 2 (L) < ■■■ < A"_(L) < -1 < 1 < A + + (L) < ■■■ < A+(L) < A+(L) 

denote the eigenvalues of L. Note that the integers M± are finite for the Jacobi matrices that appear in this 
text, but they can be zero. Similarly, let 

A] (Lk) < A 2 (Lk) < ■ ■ ■ < A k (Lr) and A^ (Lk) > A^" (Lk) > ■ ■ ■ > A k(Lk) 
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denote the real simple eigenvalues of any K x K truncation Lk, labeled in increasing (A • ) and decreasing 
order (A+), respectively. Then 


A j (L) < Xj (L k ) for j = 1,2,..., min {M , K}, 
Xj~( L) > A+ (. L k ) for j = 1,2,..., min{M + , }, 


and 

— 1 < A J{Lk) for min{M _ , K} < j < K, 

1 > AA(Lx) for min{ M + ,K} < j < K. 

In other words, pure point spectrum of L shrinks around the a.c.-spectrum under truncations. In particular, 
an eigenvalue of L that is close to the a.c.-spectrum might go inside the interval a ac (L) = [—1,1] after applying 
a truncation, in which case it fails to be captured by conventional eigenvalue algorithms that are run on the 
finite truncations Lk- The example to be discussed below illustrates such a case. For a more detailed account 
on the spectra of finite truncations of doubly-infinite Jacobi matrices, we refer the reader to Section 4 in [5] 
and . 

Consider the data {o°, 6°}: 


(65) 


b° 

u n 

f 


^ 1 y/^n—l&n+l 

~2 4 ’ 

1(1 _ 10- 4 )(z* - 1/z*) 

1 + e- 4 "/ 5 , z* = e -2 / 5 . 


4-i-n 
4-i J 


which is created by inverting pure soliton initial data. Define the ( 2K + 1) x (2 K + 1) truncation 


L2K+1 


(b-K 

a~K 

0 


\ 

a~K 

b-K- 1-1 

a~K +1 



0 

a~K +1 







a-K -2 

0 



0 >K -2 

bx-i 

ok- 1 

V 


0 

a-K-i 

bi< / 


for I\ = 2400. Using standard eigenvalue algorithms for tridiagonal symmetric matrices on L 2 K +1 yields no 
eigenvalues outside [—1,1]. However, the transmission coefficient has a pole outside [—1,1], which can be 
captured by Newton iteration to find zeros of 1 /R(z) and we find R(z*) = 0 for z* ~ 0.99982297716 which 
corresponds to an eigenvalue A ~ 1.00000001567 of L. Naturally, such an eigenvalue is difficult to capture via 
truncations because it does not emerge until K is very large. See Figure [B] for a illustration of the zero. 


Appendix C. Toda (/-function 


In this section we explicitly define the (/-function computed in Section 4.3 We use this form to determine 
the exact asymptotic form of the collisionless shock region and then discuss its computation. 


C.l. Construction of the Toda ^-function. Let Ao < 0 and po > 0 denote the real and imaginary parts of 
the stationary point, z 0 , respectively. For a, /3 G T in the upper half-plane with — 1 < Re a < Ao < Re/3 < 1, 
define 

F(a,/3)= f \V(p- a )iP- « -1 ) (P ~ P)(P~ P~ l ) + d P■ 

Jfi P 

Here the square root is defined with branch cuts on [f3, a] arc and [a -1 ,/3 _1 ] arc and asymptotics 

V(P~ a )(p- « _1 ) (P ~P)(P~ /4 1 ) ~ p 2 , p^-oo. 

Lemma C.l. Under the additional restriction Rea + Re/3 = 2Ao, F(a,/3) 6 [0,4]. Moreover, for fixed zo, F 
is a monotone decreasing function of Re a as Rea increases from —1 to Aq. 
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Figure 29. The reflection coefficient R(z) for the initial data (65). A zero for 1 /R{z) is 
observed near z = 1 corresponding to an eigenvalue of L. The zero is very close to the unit 
circle and therefore the corresponding eigenvalue is close to the right edge of the continuous 
spectrum of L. This eigenvalue is difficult to detect with traditional eigenvalue techniques but 
it is easily captured using the reflection coefficient. 


Proof. For notational simplicity set \ a = Rea, u a = arga, and \p = Re/3, up = arg/3. Also, set 

X(p) = y/{p - a) (p - a' 1 ) (p - /3) (p - /3 _1 ). 

Then 

X + (p) dp 
Ip P p' 

First, note that 

/ Y+(r h\ 2 

( 66 ) 


F(aJ3)= f 
JB 


(^7^) =^ z ~ X o t ){z-\p), 


for z = (p + p~ 1 )/ 2. Now, for p on the unit circle, z = Rep and the right hand side of (66) is negative in the 
domain of integration, vanishing only at z = \ a and 2 = A p. This implies that X + (p)/p is purely imaginary. 

Second, by definition X{p)/p ~ p as p — > oo so that Imp > 0 in the upper half-plane for \p\ large enough. 
We claim that Im(X(p)/p) cannot change sign outside the unit circle off the real axis. Suppose that it did. 
Then X(p*)/p* = c would be real valued for some p* with |p*| > 1. For z * = (p* + pf 1 )/2 

(67) 4 (z* — A a ){z* — A p) = c 2 E M. 


Solving (67) for z*, we obtain 


z ,* = 


A q + A /3 ± \J (A q — A p) 2 + c 2 


which is clearly real valued. Thus p* is on the real axis or the unit circle. This implies that X~(p)/p has a 
positive imaginary part on the arc (/3, a) arc , and hence X + {p)/p is purely imaginary with a negative imaginary 
part on (/3, a) arc . This, together with the fact that dp/p = idu in polar coordinates, p = e lu , implies that 


( 68 ) 


F(a,/3 ) = 2 / \J —(cos u — Ao ) 2 + (A a — Ao ) 2 du > 0 . 

Jup 


Setting y = cos u in ( 68 ) gives us 
(69) F(a,l 3) = 2 


f 


-1 


p V~(y - A o) 2 + (Aq, - Ao) 2 dy. 


h\ 0 -\ a \J\-y 2 

Now, for fixed Ao with A Q < Ao < 0, we differentiate F with respect to A Q : 

d L = 2(x _ X] f 2Xo ~ Xa _ ^ _ 

9K [a 0> Jx a v / i 3 7^-(y-Ao) 2 + (A a -A 0 ) 2 ' 

Here the contributions from the endpoints vanish. Since 1 — y 2 > 0 and — (y — Ao ) 2 + (A a — Ao ) 2 > 0 for 
y £ (A q ,2Ao — A q ), < 0 f° r fixed Aq. This implies that F is a monotone decreasing function of A Q as 
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A q increases from —1 to Ao- Thus, for any given —1 < Ao < 0, F attains its maximum value F*( Ao) when 
A q = —1. Moreover, 

r»2Ao+l _ ___ 

dy > 0, 


OF* 


d\r 


= 2 


i + y 


-i n/ 1 - y 2 V~(y - A o) 2 + (-i - A 0 ) 2 


which implies that F* is a monotone increasing function of Ao as Ao increases from —1 to 0. Setting zo = i, 
a = — 1, and f3 = 1 gives 

/*7T _ /*7T 

F(— 1,1)= / 2 yj 1 — (coscj) 2 du = 2 / sin cj dcz = 4, 

JO Jo 

which is the maximum value of F subject to the constraint A Q + \p = 2Ao- Clearly, the choice a = [3 = zo is 
admissible and it minimizes F at F = 0, whence we conclude F(a, j3) 6 [0,4], □ 

Similar to the case for the KdV (see 13021), introduce the variable 


s = — - 


log/Jn 


t 


Restricting to Rea + Re/3 = 2 Aq, with Ima, Im/3 > 0, the expression 


\v(p- a ) ip- a ^ (p - P) (p - 13 x ) + dp = s 


Jy r 

dehnes both a(s) and /3(s) since s is a monotone function of Rea (or, equivalently, of arga). We define the 
o-function to be 
(70) 


si?)=\ r 

2 Jp{s) 


y/[p-a)[p-a I )(p-P)(p-P x ) , 1 f a(s) y/{p-a)(p-a x ) (p - /3) {p - /3 x ) 


p z 


dp+ 2 


1-1 


p- 


dPi 


We choose the branch cut for yj (z — a)(z — /3) (z — a -1 ) ( z — f3~ l ) to contain the circular arcs on the unit 
circle from (3 to a and a~ x to (3. In order for g(z) to be single-valued, it is necessary to add a branch cut on 
the arc connecting a and a -1 . 


Lemma C.2. The g-function given by (70) satisfies: 

( 1 ) 


9 + {z) - g (z) = < 


( 2 ) 


[ \v(p~ a )(p- a X ) (p - P) (p - P l ) + dp, z£(/3(s),a(4 

Jfts) P 


[ \\J{p - a) (p - a l ){p-f3){p-fi l ) + dp, z € (a(s), a(s) x ) 

Jfy) P ^ 

f 4 ViP ~ «) iP ~ « _1 ) Ip -P)(P~ /3 _1 ) + dp, z e (a(s)~ 1 ,/3- 1 (s)) i 
■J fi(s)- 1 P 

0, 


otherwise. 


(71) g + (z)+g (z) = < 


[ \y/(p-a)(p- a” 1 ) {p - fi) (p - /3” 1 ) dp , z G (/3(s), a(s)) arc 

J-i P 

/ Z ^ 

^V(P~ a )(P~ a ~ l ) (P -P){P~ /3 _1 ) dp, z G (a(s),a(s) _1 ) are 

J ^ ^ v^(P - «) (P “ « _1 ) (p - /3) (P - /3” 1 ) dp, z G (a(s)^ 1 ,/3^ 1 (s)) ( 


2p(z), 


otherwise. 


(3) 


g(z) = -z — Aq log z + 0 (z x ) , as z —>• oo. 


(4) <?(z) — ;j70(z) bounded. 























50 


D. BILMAN AND T. TROGDON 


(5) 


ple^ 9+{z) 9 = 1, for z € (a(s),a(s) 4 ) ( 


Proof. As before, set X(p) = sj\p — a) (p — a x ) (p — /?) (p — /3 x ), and let X a and Xp denote the real parts 
of a and (3, respectively. The properties (1) and (2) follow from contour integration and the fact that 
X+ (p -1 ) = X+ p 2 P ' > . To prove (3), we use 

V 1 -y = 1 - \y - \v 2 + ° (y 3 ) » as y -> o, 

and the condition that X a + Xp = 2Ao, to obtain 

X(p) 


P 


— P 2 — 2Aq p 1 + (l — l(X a — A^) 2 ) + 4AqA q A^p + O (p 2 ) , as p —>• 0 . 


Also, for large p, we have 
X(p) 




= l-2A 0 p 1 + (1 - |(A q - Xp) 2 )p 2 + 4A 0 A q A pp 3 + O (p 4 ) , asp^oo. 


Setting H(z) = 
(72) 


l(s) P 


X(p) , . . 

— 2 — dp, integration gives us 


H(z) = —z 1 + C — 2Aologz + (l — |(A a — A^) 2 ) z + O (z 2 ) , as z —> 0, 

H(z) = z — 2Aologz + C — (1 — |(Aq, — A^) 2 ) z~ 1 + O (z~ 2 ) , as z —> oo, 

for some complex constant C. Then lim H(z) + H (^ _1 ) = 2 C, which implies, for large \z\, that 


2 C = 


(73) 


r 

J a(s) P Ja 


« P 


X(p ) , 

~^dp = 


*(s) P 


^ dp+ j^ ) -'^ dp 


P 


r “ w ^ dp = - 2 /_“ w ^) dp . 


J cn(s) 1 P 

Therefore, g(z) = \z — Aologz + O (z~ 3 ) as z — > oo. (4) follows from the integral representation and the 
asymptotic expansion around z = 0 in ([72]). Finally, to prove (5), assume z € (ct(s),a(s) _1 ) , 


.*** s _ r M x<*> ip 

t Jp{s) P 
- log pI = t (g + (z) - g-{z)) 

e t{g+(z)-g-(z)) p 2 = L 


□ 


C.2. Derivation of the collisionless shock scaling. We have 


P o = 



As described in Appendix C.l 


we choose A = K(a) and B = K(f3) (see 


for — ~ 1. 

t 

(|35[) for the transformation K) by 


(74) 


log pi . r A \J(q- b ) {q + a) (q + B) 

tpo Jb (Ao - ipoq ) 2 


Note that as Ao —> —1, po —3 > 0 (i.e. zq —> —1), the integral on the right hand side converges provided that A 
and B converge to finite values. To ensure this, we enforce that 


log Po 

Pit 


-*c. 


as t —> oo, 
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for some constant CgM. Let n = t — , with 


(75) 

Then 


0 as f —> oo, so that 
Po = '/pit) (l + 0 . 


logPo 


-log 


P(t) 

t 


pit t 

For this limit to exist, set g{t) = Cit 1 ' /3 (logt) 2,/3 , Ci > 0. Doing so yields 


(<f ) 3/2 


C , as t 


oo , 


-log ( 




) _ I log t - I log(log t) - log Cl 


(^) 3/2 


t 


C 3 / 2 logt 


3 C 


3/2 ’ 


as t oo , 


as desired: In this limit, A and B in (74), tend to finite values. Therefore the scaling for the collisionless 
shock region is given by 

n = f-Cif 1/3 (logt) 2/3 . 

C.3. Computing the Toda ^-function. To compute g(z ) we first compute g and use the relation tg(z) = 
tg(z) + ^6 (z). It follows that g^z) solves the following RH problem 

(76) g'(z) + q'(z) = t~ 1 d'(z), z £T, u U 

(77) g'( 2 :) = 0(z~ 2 ), as z —> oo. 


Furthermore, g'(^) is a bounded function on C \ (£ u US;). We remark that imposing that (76) and that g 7 (^) 
is a bounded function in the finite plane uniquely determines g'(z) and (77) is a consequence of our choice of 
a and f3. 

We consider the function &(k) = g(M(fc)), bounded in the finite plane, where M(k) = /P maps the real 
axis to the unit circle. Then &'(k) solves 

(78) 0'(A:) + 0'(jfe) = t“ 1 0'(M(A:))M , (A:), z € (- B , -A) U (A, B ), 


(79) 


= 0((k — i ) 2 ), as k —> i. 


Here A = M _1 (o:) and B = M~ 1 (/3). As in the case of g’(z) the behavior at k = i does not need to be 
imposed — the boundedness of & along with the jump condition (78) are enough to uniquely determine the 
function. Thus the numerical methodology in [26j applies directly to this situation and allows us to compute 

both ©' and © to within machine precision, uniformly in the complex plane. Hence g(z) = © (i/ij'j can be 
computed accurately. 


Appendix D. The vanishing lemma and the unique solvability of RH problems 

Let G*(z ) denote the Hermitian transpose of the matrix G(z). In the following lemma, we allow a solution 
of a RH problem to fail to be continuous up the boundary but it must be uniformly bounded and satisfy the 
jump condition almost everywhere (a.e.). 

Lemma D.l (Vanishing lemma). Consider the RH problem 

4> + (^) = <J>_(z)G( 2 ), a.e. z£ T, <h(oo) = 0, 

where is uniformly bounded on C \ T. Assume G £ L°°(T) and G*(z ) + G{z) is positive semi-definite a.e. 
and strictly positive definite on a set of positive measure. Then = 0. 

Proof. Because 4)(z) vanishes at infinity, <f>*(l /z)z~ l is a bounded analytic function on the unit disc. Then 
for \z\ = 1, 1/7 = z so that 

0 = ~i [ <5> + {z)<5>*_(l/z)z- l dz= f <S> + (z)$*_(z)\dz\ = [ <5> + (z)G*(z)$* + (z)\dz\. 

J T J T J T 
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Adding this equation to its Hermitian conjugate we find 


$+(*) [G*(z) + G{z)\ $* + (z)\dz\ = 0. 


And if G* ( z ) + G(z) is strictly positive definite on a set S then <3> + = 0 a.e. on S. As S has positive measure, 
it follows from classical results (see, for example, nu) that <f> + = 0 and therefore <f> = 0 . □ 


D.l. Unique solvability without solitons. 

Proposition D.2. The associated matrix RH problem to RH Problem [7] (i.e dropping the symmetry condition 
and normalize to I at oo) is uniquely solvable in the absence of residue condition ( i.e. if N = 0J. 

Proof. This associated matrix RH problem is uniquely solvable if the operator defined by u i-a u — CjU - (J — I) 
is invertible on H 1 ^ T) [40j. From classical results, this operator is Fredholm and from det J(z) = 1 the index 
is zero. Then Lemma |D. 1 1 demonstrates that the kernel must be trivial and the operator is invertible because 
R (.U 1 ) = R(z) and hence J(z) + J*(z) is diagonal with non-negative diagonal entries. □ 


We now introduce a function A s (z) = A s (z;zq) that satisfies RH Problem |4| with the symmetry condition 


(80) 


A{z) = (j J) (J J), 


\z\ > 1, 


with A s (oo) = C and generally C 7 ^ I. A s (z) is called the partial transmission coefficient in | 21 ] (when N = 0) 
and is given by 


A s (z) = diag (<5 s (z),<5 s H~)) , 

M z ) = ex P ( 77 — [ logr(s) 

\ 2m 


s + z ds 
s — z 2s 


, t{z) = 1 — R(z)R(z ) 


We now verify that A s satisfies the jump condition and singularity conditions of RH Problem [4] along with 
the symmetry condition (80). We begin with the jump condition. For 2 e ( 20 , 2 / — 1 using the 

Sokhotski-Plemelj lemma [231 p. 42] 

log(<5 s + (z)/8~(z)) = log t(z). 


For | 2 | > 1, 




Then sending s i—>• s 1 , using that r(s 1 )=r( 5 )we have 


5 s {l/z) = 1 /5 s (z), 


showing that A s satisfies (80). We then find 


= (irj 1 ) 


s + z 


2s 


- 1 


ds 


s — 2 


= exp 


2iri 


[ log t(s)- 


Here <$(z) = A 11 ( 2 ) where A is the solution of RH Problem |4j So S s (z) = cS(z ) for some cGC, c/0. 

Consider a new RH problem constructed using the same jump matrix J from RH Problem [2j Assume 
N = 0 (without solitons) and let C ZQ = {z : \z — zq\ = e} with counter-clockwise orientation. Let wo and vo 
be the intersection points of [zq, an d C ZQ and T \ (zq, Zq 1 )^ and C Zo , respectively. Define the jump 
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matrix J : T -> <C 2x2 , T = (w 0 ,w 0 ) arc U C Zo U C-i U (T \ (u 0 , v Q ) arc ): 


J{z) = < 


A s (z)J(z)A+(z) \ z G (wqiWq 1 ) 
A s(z)J(z), 

As(.z)j 

A s (z)J(z) A s (z) -1 


z G U (7 2 -i and |z| > 1, 
z G (7 20 U (7 -i and \z\ < 1, 


2GT \ 


It follows that J satisfies the product condition [3D] Def. 2.55] and therefore u i—>- u — C^ u- (J — I) is Fredholm 
on Hi (F) (see [4(11 . Def. 2.48]) and det J = 1 implies the index is zero. Assume u is in the kernel of this 
operator, so that 'F(z) := C^u is a solution that is continuous up to T that satisfies ’F(oo) = 0. Define for 
z G C \ (Tuf) 


'l(z) = T(z) < 


A, 


V(*). 

J~\z)A-\z), 

J-\z)A~\z), 

L 


|z — zo\ < e and \z\ < 1, 

| z — Zq 1 ] < e and \z\ < 1, 
|z — zq\ < e and \z\ > 1, 
\z — Zn 1 ] < e and \z\ > 1, 


otherwise. 

It follows that T has a continuation that is analytic in C \ T with the jump 

£+(*) = if-(z)G(z), 

'[1 - R(z)R(z- 1 )]S~(z)/S+(z) -WK(zW(z)e-<>^' 


G(z) '.= A-{z)J{z)At{z)- 1 = 


-ftp) r 9(z;n,t) 


4 + (*)M> ( z ) 


\ S a (z)5f(z) 

This follows because 'L + (z) = T _ (z) on C Zo and (7,-1. Note that by the singularity conditions in RH 


Problem Hj G(z) is continuous near z = —1. It follows that S s (z) also satisfies 6t(z) = 1 /S s (z) so that 


G*(z) = 


'[l-R(z)R(z- l )]\5~(z)\ 2 R{z)5-{z)5f{z)e- e ^)' 


_ R( z ) r 0(z;n,t) 


and 


Applying Lemma 


<5 S (z)S^(z) 

G-(z)+G(z) = (4 1 - SUWOllVMI 


I4 + (^)P 


0 

2 | 5+(^)| 2 


D.l 


'F = 0 and hence T = 0. This gives the following proposition: 


Proposition D.3. The matrix RH problem with jump matrix J is uniquely solvable. 

D.2. The addition of solitons. Our main approach to adding solitons is to just include the jumps on the 
contours in RH Problem |3j Because the matrix RH problem without these contours is uniquely solvable 
it is reasonable to expect that the addition of these jumps will not completely destroy the unique solvability 
of the problem. As stated in Remark 6.1, for each fixed n, the problem is either solvable for no t or solvable 
on the compliment of a discrete set of t values. One should expect the latter and this is indeed the case: The 
matrix RH problems we consider are uniquely solvable for sufficiently large t. 

We now present an approach that incorporate the symmetry condition (see (|11[)) for the vector RH problems 
to give equations that are uniquely solvable for every n and t value. We emphasize that this approach is not 
necessary to compute the solution and the fact that the numerical method presented here is robust despite 
ignoring symmetry is important, as evidenced in Figure [23} Ignoring symmetry is also often more convenient 
for implementation as well. Nonetheless, inspired by 121, (4.7)] and by the deformation for in RH 
Problem [3j we define 

£ _ £ . 

Qs(z) = diag(g s (z),l/g s (z)), q s (z) = ~ J 


j&K n ,i 


zQ - 1 
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We use this here because 


Q S (*) = (J JjQsCOu l 


and if m(z) (from RH Problem [lj satisfies the symmetry condition (11) then so does m*(z) := m(z)Q s (z). 
Define 


h+(n,t):=h+ = ! limz ^ Z Ci)2(?s(z) 2 ’ { G KrhU hj(n,t):=hj = \ limz ^ l{Z C ^ )2<ls ^ 2 ’ 3^ K n,t 


Qs{Cj) ! j 0 &n,t 


3 J [qsiCj)- 2 , J?K n ,u 3 

It follows that m*(z) satisfies the following conditions. 

RH Problem 10. Find the function m* : C \ T —> C lx2 that is sectionally meromorphic, continuous up to 
T, with simple poles at C^ 1 , j = 1 • • •, IV, and satisfies: 

• the jump condition: 

m+(z;n,t) = 7n k _(z-,n,t)Q~ 1 (z)J(z-,n,t)Q s (z), z£l, 

• the residue conditions: For j = 1, 2,..., N, 

Res m*(z', n, t) = lim m*(z m , n, t) I , j € K nt , 

*=Cj *~>Cj V U U J 

Res m*(z;n,t) = lim m*(z\n,t ) ( . , _i -eu r ,n,t) n ) > G K n,u 

' q 1 j *-*, -1 V ^"j, j e "/ 

Res m*(z; n, t) = lim m*{z;n,t ) ( r , n ,t)/ h + q) ’ £ Kn ^ 


Z=(j 


2->-Cj 


/n cr 1 'v- P o(Cj-, n ,t) /h~ 
Res m*(z;n,t)= lim m*(z-,n,t) ( 3 3 n 3 

*= c / 1 ^C " 1 0 


J £ Rn,t) 


the symmetry condition: 


m* (z x ’,n,t) =m*(z-,n,t) ^ ^ , 


the normalization condition: 


lim m*(z; n. t) = (to* to^) , with = 1 and (j I to* > 0 . 


Now, let K(z) = K(z\n,t) be the solution of the following RH problem: 

RH Problem 11. Find the function K : C \ T —>• C 2x2 that is sectionally analytic, continuous up to T, and 
satisfies 

• the jump condition: 

K+(z ) = K_(z)Q~ 1 (z)J(z;n,t)Q s (z ), z6l, 

• the symmetry condition: 


( 81 ) 


k (*■')=(; j) 


f/ie normalization condition: 


I\{ oo) = I. 


RH Problem 11 has a unique solution. Indeed, because q s (z) = l/q s (z) for z £ T, the jump matrix and 
contour satisfy the hypotheses of the vanishing lemma, Lemma |D.1[ This implies that the solution without 
the symmetry condition exists and is unique. Then from the symmetries of the jump matrix 

*«-(! 
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is also a solution and therefore (81) follows. Furthermore, det K(z) = 1. Define Ej t ± by 


Ej,+ < 


K (Cj 


'0 h t <, js-i 


0 


3 G K n j, 


0 


K «M_^ je e { ^ )/h + o K ^ 


Ej,- := 


_ _ *«T‘> (_ ofc - 7 -v»(c^,» o) ' 6 A ”’*, 


*(C 


-ij0 


' ' 0 0 


J /I-HC 1 ), 3 <? K n . t . 


Now consider v*(z) := m*(z)K 1 (z) which has an analytic continuation across T. Thus v*(z) satisfies the 
following discrete RH problem: 


RH Problem 12. Find the function v*: C —> C lx2 that is sectionally meromorphic with simple poles at C^ 1 , 
j = 1N, and satisfies: 

• the residue conditions: For j = 1, 2,..., N, 

Resu*(z) = lim v*(z)E j>+ , j G K n>t , 

Z=Cj z^Cj 

Res v*(z) = lim v*(z)Ej_, j G K n>t , 

^c- 1 

Resw*(z) = lim v*(z)E ji+ , j g K Ujt , 

Z=Cj z^-Cj 

Res v*(z) = lim v*(z)Ej , j 0 K nt , 


2 =C7 




• the symmetry condition: 

(82) u* (z- 1 ) = v*{z) (J j), 

• the normalization condition: 


(83) 


lim v*(z) = (u* v%) , with v* ■ = 1 and J Q J v\ > 0 . 


We assume we can solve for K(z) and therefore compute each of Ej t ±, j = 1, 2,..., N. It follows that 

= (c + Ef.i«,.+$£ + Ef.i d + Ef.,%+S^ + Ef., h- 0 h ). 


for some choice of constants a J: ±. b^±. From the symmetry condition (82) it follows that aj.± = bj - F , c = d 
so that 


”*M = (c + EjEj + Ef., c+ Ef., <*-*£ + Ef=, %,+ igh, . 

We now obtain a linear system for these constants under the assumption that c is known. For k = 1,2,... ,N 


lies v*(z) = [a k ,+ {(k ~ !) «fc,-(Cfc - !)) 

—Cfc 

(84) = (c + E^%,+^ + Ef.,%-^ c + E^%-|^ + Ef. 1 % + fcSr)^ + . 

Res u*(z) = (a fe -(Cfc 2 - 1) afc,+ (C 2 - 1)) 

2 =C*, 




Cfc'O-i 


cr 1 -^- 


if? 


CF 1 ^-! 


1 + aj ’ + CF i C 3 --i ) Efe -- 


(85) 
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This is a system of 4 N equations for 2N unknowns so it must have redundancy. We compute the symmetries 
of Ej ± from the symmetry of v*(z): 


Resu*(z) = lim v*(z)Ej + , 
z=Cj z->£j 

lim (z — (j)v*(z) = lim v*(z)Ej j +, 



This leads to the conclusion that (84) is equivalent to (85). And since a unique solution of (84) exists for 
the correct c, a solution must exist for every c. Due to linearity, the coefficients a,j t ± = (ij,±(c) have simple 
dependence on c. Namely, a h ±{c) = ^±(1) • c. To find c, we set c = 1 and solve for a,j t ±( 1). Then for 
v*(z) = v*(z;c ) we have 


And c is then found by enforcing (83). 


u*(0; c) = cu*(0;l). 


Appendix E. A proof that generically R(— 1) = -1 

In this appendix we prove a theorem to establish the genericity of i?(±l) = — 1 for the Jacobi matrices used 
in this work. Consider the weighted £ 1 -space of doubly infinite sequences, £^(Z), with the weight function 
given by n i-A 1 + |n|, and define the Banach space £^(Z) = ^(Z) © £^(Z) equipped with the norm 

\\i.x,y)\\ci = + M)(M + \yn\) ■ 

nEZ 

We let A4 denote the Marchenko class of doubly infinite Jacobi matrices J(a, b) whose coefficients have the 
property that 

({ a n — 2 Jnez ’ {bn}neZ^ ^-E W (Z). 

For ease of notation, we set a n = a n — \ for each n E Z and define J 0 (-, •) by J Q (a — 1/2, b ) := J(a, b) for all 
of the Jacobi matrices mentioned throughout this appendix. Note that J 0 (a,b) = J(a,b). 

Theorem 1. The set of doubly infinite Jacobi matrices J 0 {a, b ) in A4 with the associated reflection coefficient 
satisfying R(± 1) = —1 is an open and dense subset of Ai in the topology induced by the norm 

\\Jo(a,b)\\ M = X ( 1 + N) (W + IM)- 

nEZ 

Proof. We first show that the subset of Jacobi matrices in A4 with 12(4=1) = —1 is dense in A4. Suppose that 
L = J 0 (a, b ) is in M with the reflection coefficient satisfying Rl{- tl) ^ —1, and let e > 0 be given. Also, let 
Lo denote the free Jacobi matrix with coefficients a n = \ ( i.e. a n = 0) and b n = 0 for all n E Z - namely, 
the discrete Schrodinger operator with the zero potential. As J 0 (a, b) £ A4, there exists some N = N(s) £ N 
such that 

X (! + M)(M + \bn\) < 

\n\>N 


(86) 
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Now consider L = J 0 (a,b) which is the Jacobi matrix that is a finite-rank perturbation of Lq, with the 
property that 




a n , if — N < n < N — 1, 
0, otherwise, 

Define the Jost solutions ip± and p\ t of the problem 


and b n = 


if — N < n < N, 


0, otherwise. 


(87) 

by their asymptotic behaviors: 

( 88 ) 


v -i 


~ z H - z 

( Lp){n ) =--- p(n), for all n € Z, 


lim p\ (z;n)z n = 1 and lim p T _(z-,n)z n = 1, 

n—>+oo n—>-+oo 

lim p\(z; n)z~ n = 1 and lim p l _(z\n)z n = 1. 


First recall that the Wronskian of any two solutions p and ip to the problem (87) is given by 
W (p(z; -),ip(z; •)) = a n [p(z\ n)ip(z;n + 1) - p{z\ n + 1 )ip(z]n)\ = - det 


p(z;n + 1) ip(z;n + 1) 
a n p(z\ n) a n ip(z;n) 


and that it is independent of n as long as / and g solve (87) with the same value of \ (z + z _1 ). Recall also 
that the reflection coefficient Ri(z) associated with L (or for any Jacobi matrix L in Af) can be obtained by 
ratio of two Wronskians of certain Jost solutions that correspond to the same value of ^ {z + z _1 ): 

W (pL(z; -),p r (z; ■)) 


(89) 


R-l(.*) = 


W ( p r + (z;-),p 1 _(z ;■)) 


For any solution p of (87), observe that 


p(z;n + 1) 
a n p(z;n) 

= A n (z) 

p(z;n ) 

a n -ip{z;n- 1) 

and 

<p(z;n) 

a n -ip(z;n- 1) 

= A n (z) 1 

p(z;n + 1) 
a n p(z\n ) _ 


where A n (z) =tl are the transfer matrices defined by 


A-n{z) — 


^ — h — 1 

2 u n L 

at 0 


and A n (z) 1 = — 


0 


1 


—a; 


2 z+z~ 


-W, 


Now, since n i-a z n and n *—> z~ n solve (87) with L replaced with the free matrix Lq, as long as the Jost 
solutions of (87) are outside the support of L — Lq, they coincide with the solutions of the problem for the 
free matrix Lq. In particular: 

<P+(z;N + 1) = z N+1 , p+(z;N) =z N , p T _(z;N+ 1) = z~ N ~ 1 , p r _(z] N) = z~ N , 

p l + (z;—N — l) = z~ N ~ 1 , p\{z\ —N) = z ~ N , p l -(z-, — N — 1) = z N+1 , p l _{z\—N) = z N . 

To compute the reflection coefficient, we evaluate the Wronskians in (|89|) at n = 0. Then we have 


(91) 


Rd*) = 


W(p l _(zp 

),P-{z;-)) 


pl(z-,l) 

p r -(z;l) 


a o p l (z;0 ) 

aop T {z\ 0) 

W (p+(z', • 

), P-(z-, ■)) 

det 

P+(z] 1) 

aop'p^O) 

p X -(z\l) 
aop l _{z\ 0)_ 


Using the transfer matrices, this could be expressed as: 

w ' 


det 


Rl(z) = 


A 0 (z)A_i(z) • • • A-n{z) 


12 


z 

1 UV+1 


.-(JV+1)' 


1 ~—N 


l 2' 


det 


Ai(z) x A 2 {z) 1 ---A n (z) 


\-i 


z n+y 

1 ~N 
L 2* . 


; A 0 (z)A_i(z) • • • A - N (z) 


z N ' 
l r AT+l 
L2~ 


For ease of notation, label the products of transfer matrices that appear above: 

A n (z) := A 0 (z)A_i( 2 :) • • • A - N (z), 

Biv(z) := Ai(^) _1 A 2 (z) _i • • • Atv(^)' 1 , 
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and rewrite 


det 


R~M = 


A n (z) 


z N ' 
±?N+1 
L 2 “ 


; B n(z) 


~ z -(N+iy 

ly~N 
L 2 Z 


det 




n+y 


1 


12 Z 


A A r(z) 


Z N ' 
1 ~TV+1 
L2~' 


A direct calculation shows that 


holds unless 


R;(±l) = -1 


det 


Bjv(±1) 


Ajv(±1) 


1 

l 

L 2 J J 


= 0. 


In case the determinant in (|Ej) is nonzero, since (86) implies that ||L — L\\j^ < |, we have proven that the 
subset of Jacobi matrices in M whose reflection coefficients satisfies i?(±l) = —1 is dense in A4. The case 
where the determinant in (Je]) vanishes requires more work and will be treated below. For notational brevity, 
we present the argument for z = —1 only. The argument for the case z = 1 is identical. 

Suppose that (Je]) holds, and note that the determinant on the left hand side of ([E]) is a rational function 
in 2N + 1 real variables: o_tv, ..., ajy-i, ft- at, ..., bjy- More precisely, it is of the form: 


det 


Bat(-I) 


1 

l 

L2J 


; Ajv(-l) 


l 

i 

L 2 J 


y ( CL-N, ■ ■ ■ , ajv-l, ft- 


N, 


> ftAf) 


a~N ■ 


1 «o 


■ a-N-i 


where the denominator is nonzero, and y is a polynomial in 2N + 1 real variables. Therefore determinant 
vanishes at a point in M 2JV+1 if and only if y vanishes at that point. Note that y is not the zero polynomial since 
y(0,0,..., 0) = Now suppose that y vanishes at a point (a_Ar,..., a^-i, b-N, ■ ■ ■, ftiv)- Any polynomial 
that vanishes on an open set is identically zero, and since y is not identically zero, there must be a point where 
y is non-zero in every neighborhood of where it vanishes. In particular, there is a point in any neighborhood 
of (a_ 7 v,..., cln-i, ft- n, • • •, 6jv) where y does not vanish. Take the open ball D e ^ (in M 2JV+1 ) centered at 
(a - N ,..., oat-i, b - N ,..., b N ) with radius 6(1+jV) £ (1+2JV) ■ There exists a point (a '_ N ,..., a' N _ { .b'_ Nl ... ,b' N ) e 
D £j n where 

U ( a -AD • • • ! a AT-l) ft—TV> • • • j ft]v) / 0- 

This immediately implies that the reflection coefficient Rl* for the Jacobi matrix L* = J 0 (a*,b*) with the 
coefficients 

if - N <n< N - 1, 
otherwise, 


a r 


0, 


and ft* = 


bL 


0, 


if — N < n < N, 
otherwise. 


satisfies 

Moreover, 

L* -L 


R l * (—1) = —l- 


M 


+ l n D 


An — O'* I + 


nGZ 


bn. - bt 


< 


6(1 + N)(2N + 1) 


2(1 + N)(2N + 1) = -. 


It follows from (86) and (M that ||L — L*||y^ < e. This shows that the subset of Jacobi matrices in M whose 


reflection coefficients have the property that i?(±l) = —1 is dense in M. 

We now prove that such matrices form an open set in A4. When z = ±1 we note that <p± = and when 
considering ( [89] ) we see that J2(±l) / —1 only when W(tp L(±l, -),i^L(il, •)) = 0- So, we take L = J 0 (a,ft) 
such that W(pL(±l, •); ■)) 7^ 0 an d if is enough to show that this Wronskian is continuous with respect 

to the topology of M. We will show that the following mappings 


T r : M —7 £°°(Z, + ), T r (L) = cpl(±I,-), 

Tl :M^ £°°(Z~ U {!}), T,(L) = ^_(±1, ■), 
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are well-defined and continuous which imply the continuity of the Wronskian at n = 0. It follows that 

OO 

</_(±l,n) = (±l) n - Y, (±1 ) m - n (m-n)/4, 


m=n -\-1 
71—1 




(±l,n) = (±l)”- ]T (±l)’”-"("> - «)C, 


771= —OO 


h T m = a m -i<P- (±1, m — 1) + b m ip T _(± l,m) + a m <+(±l,m + 1), 
h l m = 1) + m) + a m <+(±l, m + 1). 


+_(±l,n) + K T n m {a, b)(p r _(±l,m) = (±l) n , 

771=71+1 

71—1 

<+(±l,n) + Y K n,m( a > b )'P 1 -(± 1 >m) = (±l) n . 


We do not construct the kernels in the sums explicitly, but rather, obtain bounds on them. In both cases 
Ini < |m| so that 


We then write 
(92) 


sup Y K n,m( a , b ) < 8\\ J o(a,b)\\ M , 

n>0 m=n+1 

71—1 

sup Y K li,m( a > b ) < 8\\ J o(a,b)\\ M - 

n < 1 m=—oo 

Then it follows from [531 Lemma 7.8] that <+ and tp[_ are the unique solutions of these Volterra summation 
equations. Now, let / be the solution of (92) with a different Jacobi matrix J 0 (d, b). Then u := ip r _ — f solves 

OO OO 

u{n) + Y K n,m{a,b)u(m) = Y ( K n,m{ a ~ «» b ~ &))+_(±l, 17l). 

771=71+1 771=71+1 

Then by [Ml Lemma 7.8] 

IMI*°°(Z+) < 8||^-II^(Z+) \\Jo(a — a,b — 6)1+exp (8|| J 0 (d,6)|+) 

and hence if J 0 (a,b ) —> J Q (a,b ) in the topology of M then / —> ip T _ uniformly. Similar arguments follow for 
the continuity of <+ in this topology. This completes the proof. □ 

E.l. A simple example. Suppose that b n = 0 for all n ^ 0, 6o = /3, and we have a n = ^ for all n E Z, then 
y4(z;l) = «, <+(z;0) = l, <f?_(z;l) = z~ 1 , <+ + 0) = l 


(93) 


Using these gives 


ip\(z;l) = z-2/3, ip l + (z; 0) = 1, <+ + 1) = z 1 - 2/3, <+(z;0) = l. 


det 


(94) 


i?(z) = 


+ (u i) <++i) 
L 5<+ + 0) ^L + O) 


2/3 


det 


¥>+(*; i) +(-;i) 

}<p+(z; 0) |<+ + 0) 


2 ” 1 -z-2/3 


If (3 ^ 0, then i?(=tl) = —1. If j3 = 0, then we have the free matrix Lo, which has R(z) = 0 for all z on the 
unit circle. 
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